Code to reproduce results of the manuscript ‘Kidney Failure Prediction: Multicenter External Validation of the KFRE Model in Patients with CKD Stages 3-4 in Peru’

Author

Percy Soto-Becerra

Published

March 29, 2023

1 Introduction

This document presents the code and results of the main analysis shown in the article.

2 Setup

Mostrar código
rm(list = ls())

# Use pacman to check whether packages are installed, if not load
if (!require("pacman")) install.packages("pacman")
library(pacman)

# Unload all package to begin in a session with only base packages
pacman::p_unload("all")

# Install packages
pacman::p_load(
  here, 
  skimr, 
  survival,
  rms,
  cmprsk,
  riskRegression,
  mstate,
  pseudo,
  pec,
  riskRegression,
  plotrix,
  knitr,
  splines,
  kableExtra,
  flextable,
  gtsummary,
  boot,
  tidyverse,
  rsample,
  gridExtra,
  webshot, 
  patchwork,
  survival, 
  cmprsk, 
  survminer, 
  ggsci, 
  cowplot, 
  scales, 
  patchwork, 
  labelled, 
  glue, 
  dcurves, 
  broom, 
  downlit, 
  xml2, 
  gghalves, 
  devtools, 
  htmltools, 
  gghalves, 
  ggtext, 
  DiagrammeR, 
  gt, 
  janitor
)

if (!require("smplot2")) devtools::install_github('smin95/smplot2', force = TRUE)

# Import data
data <- readRDS(here::here("Data/Derived/data_derived.rds")) 

# Subset patients with CKD Stages 3a-3b-4
data %>%  
  filter(ckd_stage == "Stages 3-4") -> dataA

# Subset patients with CKD Stages 3b-4
data %>%  
  filter(ckd_stage2 == "Stages 3b-4") -> dataB

3 Descriptive analysis

In total, 7519 patients were included in the analysis of the group with stages 3a-4 CKD, and 2798 in the group with stages 3b-4. All patients had outcome data. The number of events was greater than 100 in all populations and outcomes, except for renal failure at 2 years in the population with stages 3b-4 (n = 88), so estimates in this group should be interpreted with more caution. Specifically, in the 3a-3b-4 group, 114 (1.5%) patients developed renal failure at 2 years and 239 (3.2%) at 5 years. Many more patients died without experiencing renal failure: 563 (7.5%) at 2 years and 1400 (18.6%) at 5 years.

Regarding the group restricted to stages 3b-4, 88 (3.1%) patients developed renal failure at 2 years and 182 (6.5%) at 5 years. Similarly, many more patients died without experiencing renal failure: 300 (10.7%) at 2 years and 683 (24.4%) at 5 years.

The median observation time was 4.9 years and the maximum follow-up was 7.8 years in the 3a-4 group. The median and maximum observation times were similar in the 3b-4 group, at 4.9 and 7.8 years, respectively.

3.1 Fig 1

Mostrar código
# Create grid of 100 x 100----
data_flow <- tibble(x = 1:100, y = 1:100)

data_flow  %>% 
  ggplot(aes(x, y)) + 
  scale_x_continuous(minor_breaks = seq(1, 100, 1)) + 
  scale_y_continuous(minor_breaks = seq(1, 100, 1)) + 
  theme_linedraw() -> p

# Create boxes----
# 

box_xmin <- 33 - 20
box_xmax <- 75 - 20
box_ymin <- 94
box_ymax <- 99
box_size <- 0.25

text_param <- function(box_min, box_max) {
  mean(c(box_min, box_max))
}

text_size <- 2.5


p + 
  # Col 0----
  ## Level 1----
  geom_rect(xmin = box_xmin, xmax = box_xmax, 
            ymin = box_ymin - 1, ymax = box_ymax + 1, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin, box_xmax), 
           y = text_param(box_ymin - 1, box_ymax + 1), 
           label = 'Total patients in VISARE database\n(n = 57,308)', 
           size = text_size ) + 
  ## Level 2----
  geom_rect(xmin = box_xmin + 5, xmax = box_xmax - 5, 
            ymin = box_ymin - 41, ymax = box_ymax - 39, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin + 5, box_xmax - 5), 
           y = text_param(box_ymin - 41, box_ymax - 39), 
           label = 'Total of patients with CKD\n(n = 22,744)', 
           size = text_size ) +   
  # Col -1----
  geom_rect(xmin = box_xmin - 13, xmax = box_xmax - 27, 
            ymin = box_ymin - 81, ymax = box_ymax - 79, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin - 13, box_xmax - 27), 
           y = text_param(box_ymin - 81, box_ymax - 79), 
           label = 'Patients with CKD 3a-4\n(n = 7,519)', 
           size = text_size ) + 
  # Col +1----
  ## Level 1----
  geom_rect(xmin = box_xmin + 23, xmax = box_xmax + 47, 
            ymin = box_ymin - 19, ymax = box_ymax - 12, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin + 23, box_xmax + 47), 
           y = text_param(box_ymin - 19, box_ymax - 12), 
           label = '34,564 screened patients did not CKD in any stage', 
           size = text_size )  + 
  ## Level 2----
  geom_rect(xmin = box_xmin + 23, xmax = box_xmax + 47, 
            ymin = box_ymin - 19 - 45, ymax = box_ymax - 12 - 45, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin + 23, box_xmax + 47), 
           y = text_param(box_ymin - 19 - 45, box_ymax - 12 - 45), 
           label = 'Not elegible (n = 15,225) \n 8,426 lacks of data on urine albumin and urine creatinin \n 6,799 were in other CKD stages (1, 2, or 5)', 
           size = text_size )  + 
  ## Level 3----
  geom_rect(xmin = box_xmin + 27, xmax = box_xmax + 13, 
            ymin = box_ymin - 81 , ymax = box_ymax - 79, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin + 27, box_xmax + 13), 
           y = text_param(box_ymin - 81, box_ymax - 79), 
           label = 'Patients with CKD 3b-4\n(n = 2,798)', 
           size = text_size )  + 
  # vertical arrow
  geom_segment(x = text_param(box_xmin, box_xmax), xend = text_param(box_xmin, box_xmax), 
               y = box_ymin - 1, yend = box_ymax - 39, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt",
               arrow = arrow(length = unit(1, "mm"), type = "closed")) +  
  geom_segment(x = text_param(box_xmin, box_xmax), xend = text_param(box_xmin, box_xmax), 
               y = box_ymin - 41, yend = 25, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt") + 
  # horizontal arrow 1-->
  geom_segment(x = text_param(box_xmin, box_xmax), xend = box_xmin + 23, 
               y = text_param(box_ymin - 19, box_ymax - 12), yend = text_param(box_ymin - 19, box_ymax - 12), 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt", 
               arrow = arrow(length = unit(1, "mm"), type = "closed")) + 
  # horizontal arrow 2-->
  geom_segment(x = text_param(box_xmin, box_xmax), xend = box_xmin + 23, 
               y = text_param(box_ymin - 19 - 45, box_ymax - 12 - 45), yend = text_param(box_ymin - 19 - 45, box_ymax - 12 - 45), 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt", 
               arrow = arrow(length = unit(1, "mm"), type = "closed")) + 
  # horizontal segment --
  geom_segment(x = text_param(box_xmin - 13, box_xmax - 27), xend = text_param(box_xmin + 27, box_xmax + 13), 
               y = 25, yend = 25, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt") + 
  # vertical arrow -->
  geom_segment(x = text_param(box_xmin - 13, box_xmax - 27), xend = text_param(box_xmin - 13, box_xmax - 27), 
               y = 25, yend = box_ymax - 79, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt", 
               arrow = arrow(length = unit(1, "mm"), type = "closed")) + 
  # vertical arrow -->
  geom_segment(x = text_param(box_xmin + 27, box_xmax + 13), xend = text_param(box_xmin + 27, box_xmax + 13), 
               y = 25, yend = box_ymax - 79, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt", 
               arrow = arrow(length = unit(1, "mm"), type = "closed")) + 
  theme_void() + 
  theme(plot.background = element_rect(fill = "white")) -> plot_flowchart

ggsave(filename = "plot_flowchart.png", 
       plot = plot_flowchart, 
       device = "png", 
       path = here("Figures"), 
       scale = 1, 
       width = 12, 
       height = 12, 
       units = "cm", 
       dpi = 600)
Mostrar código
knitr::include_graphics(here("Figures/plot_flowchart.png"))

3.2 Table 1

Table 1 provides a description of the study population. Table S1 (see section “tableS1”) details the characteristics of the study population with stages 3a-3b-4 according to the occurrence of the outcome of interest, renal failure at 2 and 5 years. The number of events of interest at 2 years was reported, and Table S2 (see section “tableS2”) shows the same information for the 3b-4 population. Table S3 (see section “tableS3”) describes the characteristics according to stages 3a, 3b, and 4 separately. With this stratification, the numbers of cases of renal failure at 2 years were very low for the subpopulations with stages 3a (n = 26), 3b (n = 36), and 4 (n = 52). Likewise, the numbers of cases at 5 years were very low for the subpopulations with stages 3a (n = 57), 3b (n = 81), and 4 (n = 101). Therefore, it was not reliable to perform predictive performance evaluation in these subgroups.

Mostrar código
dataA %>% 
  dplyr::select(sex, age, hta, dm, acr_cat, grf_cat, ckd_class, crea,
                eGFR_ckdepi, acr, urine_album, urine_crea,  death2y, 
                eventd2ylab, death5y, eventd5ylab) %>% 
  gtsummary::tbl_summary(
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({p25} - {p75})", 
    "{min} - {max}"
  ),
  digits = list(all_continuous() ~ c(1, 1, 1, 1), 
                all_categorical() ~ c(0, 1)), 
  missing_text = "Missing"
  ) %>% 
  # add_p() %>% 
  bold_labels() -> tab1a

dataB %>% 
  mutate(
    grf_cat = droplevels(grf_cat), 
    ckd_class = droplevels(ckd_class)
  ) %>% 
  set_variable_labels(
    grf_cat = "GFR categories", 
    ckd_class = "CKD KDIGO classification"
  ) %>% 
  dplyr::select(sex, age, hta, dm, acr_cat, grf_cat, ckd_class, crea,
                eGFR_ckdepi, acr, urine_album, urine_crea,  death2y, 
                eventd2ylab, death5y, eventd5ylab) %>% 
  gtsummary::tbl_summary(
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({p25} - {p75})", 
    "{min} - {max}"
  ),
  digits = list(all_continuous() ~ c(1, 1, 1, 1), 
                all_categorical() ~ c(0, 1)), 
  missing_text = "Missing"
  ) %>% 
  # add_p() %>% 
  bold_labels() -> tab1b

tbl_merge(list(tab1a, tab1b), 
          tab_spanner = c("**CKD Stages 3a-3b-4**", "**CKD Stages 3b-4**")) %>% 
  modify_caption("Table 1. Baseline characteristics of the study population according CKD Stages") %>% 
  bold_labels() -> tab1

tab1 %>% 
  as_flex_table() %>% 
  flextable::save_as_docx(path = here("Tables/Table1.docx"))

tab1 %>% 
  gtsummary::as_kable_extra() %>%
  kableExtra::kable_styling("striped")
Table 1. Baseline characteristics of the study population according CKD Stages
CKD Stages 3a-3b-4
CKD Stages 3b-4
Characteristic N = 7,519 N = 2,798
Sex
Female 4,107 (54.6%) 1,398 (50.0%)
Male 3,412 (45.4%) 1,400 (50.0%)
Age (years)
Mean (SD) 74.0 (10.2) 75.6 (10.6)
Median (IQR) 75.0 (68.0 - 81.0) 77.0 (70.0 - 83.0)
Range 23.0 - 97.0 23.0 - 97.0
Hypertension 4,486 (59.7%) 1,636 (58.5%)
Diabetes Mellitus 1,845 (24.5%) 674 (24.1%)
Persistent albuminuria categories
A1 4,772 (63.5%) 1,494 (53.4%)
A2 2,018 (26.8%) 860 (30.7%)
A3 729 (9.7%) 444 (15.9%)
GFR categories
G3a 4,721 (62.8%)
G3b 2,207 (29.4%) 2,207 (78.9%)
G4 591 (7.9%) 591 (21.1%)
CKD KDIGO classification
Moderately increased risk 3,278 (43.6%)
High risk 2,460 (32.7%) 1,302 (46.5%)
Very high risk 1,781 (23.7%) 1,496 (53.5%)
Serum Creatinine (mg/dL)
Mean (SD) 1.4 (0.4) 1.7 (0.5)
Median (IQR) 1.3 (1.1 - 1.5) 1.6 (1.4 - 1.9)
Range 0.8 - 3.9 1.1 - 3.9
eGFR using CKD-EPI, ml/min/1.73m2
Mean (SD) 46.2 (9.8) 35.7 (7.3)
Median (IQR) 48.7 (40.4 - 53.8) 37.3 (31.4 - 41.7)
Range 15.0 - 60.0 15.0 - 45.0
Albumin-to-creatinine ratio, mg/g
Mean (SD) 248.6 (3,044.4) 334.1 (3,050.8)
Median (IQR) 14.6 (4.5 - 66.1) 26.0 (6.5 - 153.8)
Range 0.0 - 144,870.6 0.0 - 144,870.6
Urine Albumin (mg/ml)
Mean (SD) 8.3 (28.3) 13.1 (36.3)
Median (IQR) 1.0 (0.3 - 4.0) 1.7 (0.4 - 10.1)
Range 0.0 - 658.0 0.0 - 658.0
Urine Creatinine (mg/dl)
Mean (SD) 72.4 (47.5) 71.4 (47.0)
Median (IQR) 63.3 (41.3 - 86.0) 64.9 (43.3 - 85.0)
Range 0.1 - 722.1 0.7 - 722.1
Death at 2 years 640 (8.5%) 358 (12.8%)
Outcome at 2 years
Alive w/o Kidney Failure 6,842 (91.0%) 2,410 (86.1%)
Death w/o Kidney Failure 563 (7.5%) 300 (10.7%)
Kidney Failure 114 (1.5%) 88 (3.1%)
Death at 5 years 1,539 (20.5%) 784 (28.0%)
Outcome at 5 years
Alive w/o Kidney Failure 5,880 (78.2%) 1,933 (69.1%)
Death w/o Kidney Failure 1,400 (18.6%) 683 (24.4%)
Kidney Failure 239 (3.2%) 182 (6.5%)
1 n (%)

4 Cumulative incidence function for competing risks data

Figure Figure 1 shows the cumulative incidence curves of renal failure and pre-renal failure death in the study patients.

Mostrar código
# Seleccion del grupo 3a-4----
cuminc(ftime = dataA$time5y, 
       fstatus = dataA$eventd5ylab, 
       cencode = "Alive w/o Kidney Failure") -> cif

ciplotdat <- 
  cif %>% 
  list_modify("Tests" = NULL) %>% 
  map_df(`[`, c("time", "est", "var"), .id = "id") %>% 
  mutate(id = recode(
    id, 
    "1 Death w/o Kidney Failure" = "Death w/o Kidney Failure", 
    "1 Kidney Failure" = "Kidney Failure"), 
    ll = est - 1.96 * sqrt(var), 
    ul = est + 1.96 * sqrt(var)
    ) %>% 
  rename(
    event = id
  )

ciplotdat %>% 
  ggplot(aes(x = time, y = est)) +
  geom_ribbon(aes(ymin = ll, ymax = ul, fill = event), 
              alpha = 0.25, linetype = 0) + 
  geom_step(lwd = 0.5, aes(color = event)) +
  theme_survminer() +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + 
  labs(x = "Years", 
       y = "Cumulative incidence",
       title = "A CKD Stages 3a-3b-4") + 
  theme(legend.position = "top",
        legend.title = element_blank(), 
        legend.background = element_rect(fill = "white"), 
        legend.key.size = unit(0.2, "cm")) -> g1

ciplotdat %>% 
  filter(event == "Kidney Failure") %>%
  ggplot(aes(x = time, y = est)) +
  geom_ribbon(aes(ymin = ll, ymax = ul), fill = "#89E1E3", 
              alpha = 0.1, linetype = 0) + 
  geom_step(lwd = 0.5, color = "#89E1E3") +
  theme_survminer() +
  ylim(c(0, 0.025)) +
  scale_y_continuous(labels = scales::percent, limits = c(0, 0.04)) + 
  labs(x = "", 
       y = "",
       title = "") -> g2

kf_fit <- survfit(
  Surv(time5y, ifelse(eventd5y != 0, 1, 0)) ~ 1, 
  data = dataA
)

num <- ggsurvplot(
  fit = kf_fit, 
  risk.table = "nrisk_cumevents", 
  risk.table.y.text = FALSE,
  risk.table.y.text.col = FALSE, 
  tables.y.text = FALSE, 
  tables.y.text.col = FALSE, 
  ylab = "Years",
  risk.table.fontsize = 3.2,
  tables.theme = theme_survminer(font.main = 10)
  )

cowplot::plot_grid(
  g1, 
  num$table + theme_cleantable(), 
  nrow = 2, 
  rel_heights = c(4, 1), 
  align = "v", 
  axis = "b"
  ) -> g3
  
g3 + inset_element(g2, 0.15, 0.43, 1, 0.856,  align_to = 'full',  
                   ignore_tag = TRUE) -> plot_cif_mh
 
# Seleccion del grupo 3b-4----

cuminc(ftime = dataB$time5y, 
       fstatus = dataB$eventd5ylab, 
       cencode = "Alive w/o Kidney Failure") -> cif

ciplotdat <- 
  cif %>% 
  list_modify("Tests" = NULL) %>% 
  map_df(`[`, c("time", "est", "var"), .id = "id") %>% 
  mutate(id = recode(
    id, 
    "1 Death w/o Kidney Failure" = "Death w/o Kidney Failure", 
    "1 Kidney Failure" = "Kidney Failure"), 
    ll = est - 1.96 * sqrt(var), 
    ul = est + 1.96 * sqrt(var)
    ) %>% 
  rename(
    event = id
  )

ciplotdat %>% 
  ggplot(aes(x = time, y = est)) +
  geom_ribbon(aes(ymin = ll, ymax = ul, fill = event), 
              alpha = 0.25, linetype = 0) + 
  geom_step(lwd = 0.5, aes(color = event)) +
  theme_survminer() +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) + 
  labs(x = "Years", 
       y = "Cumulative incidence",
       title = "B CKD Stages 3b-4") + 
  theme(legend.position = "top",
        legend.title = element_blank(), 
        legend.background = element_rect(fill = "white"), 
        legend.key.size = unit(0.2, "cm")) -> g1

ciplotdat %>% 
  filter(event == "Kidney Failure") %>%
  ggplot(aes(x = time, y = est)) +
  geom_ribbon(aes(ymin = ll, ymax = ul), fill = "#89E1E3", 
              alpha = 0.1, linetype = 0) + 
  geom_step(lwd = 0.5, color = "#89E1E3") +
  theme_survminer() +
  ylim(c(0, 0.15)) +
  scale_y_continuous(labels = scales::percent, limits = c(0, 0.08)) + 
  labs(x = "", 
       y = "",
       title = "") -> g2

kf_fit <- survfit(
  Surv(time5y, ifelse(eventd5y != 0, 1, 0)) ~ 1, 
  data = dataB
)

num <- ggsurvplot(
  fit = kf_fit, 
  risk.table = "nrisk_cumevents", 
  risk.table.y.text = FALSE,
  risk.table.y.text.col = FALSE, 
  tables.y.text = FALSE, 
  tables.y.text.col = FALSE, 
  ylab = "Years",
  risk.table.fontsize = 3.2,
  tables.theme = theme_survminer(font.main = 10)
  )

cowplot::plot_grid(
  g1, 
  num$table + theme_cleantable(), 
  nrow = 2, 
  rel_heights = c(4, 1), 
  align = "v", 
  axis = "b"
  ) -> g3

g3 + inset_element(g2, 0.15, 0.51, 1, 0.856,  align_to = 'full',  
                   ignore_tag = TRUE) -> plot_cif_vh

(plot_cif_mh / plot_cif_vh) + 
  plot_annotation(tag_levels = 'A') -> plot_cif

ggsave(filename = "Plot_CIF.png", 
      plot = plot_cif, 
      device = "png", 
      path = here("Figures"), 
      dpi = 300, 
      scale = 2, 
      width = 8.5,
      height = 14, 
      units = "cm", 
      bg = "white")

4.1 Fig 2

Mostrar código
plot_cif

Figure 1: Cumulative incidence curves of the observed outcome probabilities in the population study for kidney failure or deathd without kidney failure

5 Predictive Performance

5.1 Calibration

Mean calibration: OE ratio

Mostrar código
# Calibration (O/E) -------------------------------------------------------

# Seleccion del grupo: Stages 3-4----
vdata <- dataA %>% 
  select(id, risk2y, risk5y, time5y, eventd5y, time, eventd) %>%  
  drop_na() %>%  
  mutate(status = factor(eventd5y, 
                         levels = c(0, 1, 2), 
                         labels = c("Alive w/o Kidney Failure", 
                                    "Kidney Failure", 
                                    "Death w/o Kidney Failure")))
primary_event <- 1

# A 2 años----
horizon <- 2

vdata$pred <- vdata$risk2y

# First calculate Aalen-Johansen estimate (as 'observed')
obj <- summary(survfit(Surv(time5y, status) ~ 1, 
                        data = vdata), 
                times = horizon)

aj <- list(
  "obs" = obj$pstate[, primary_event + 1],
  "se" = obj$std.err[, primary_event + 1]
)


# Calculate O/E
OE <- aj$obs / mean(vdata$pred)

# For the confidence interval we use method proposed in Debray et al. (2017) doi:10.1136/bmj.i6460
k <- 2
alpha <- 0.05
OE_summary <- cbind(
  "OE" = OE,
  "Lower .95" = exp(log(OE - qnorm(1 - alpha / 2) * aj$se / aj$obs)),
  "Upper .95" = exp(log(OE + qnorm(1 - alpha / 2) * aj$se / aj$obs))
)

OE_summary2a <- round(OE_summary, k)

avg_obs <- cbind(
  "avg_obs" = aj$obs * 100,
  "Lower .95" = 100 * (aj$obs - qnorm(1 - alpha / 2) * aj$se),
  "Upper .95" = 100 * (aj$obs + qnorm(1 - alpha / 2) * aj$se)
)

avg_pred <- cbind(
  "avg_pred" = mean(vdata$pred) * 100
)

avg_obs2a <- round(avg_obs, k)
avg_pred2a <- round(avg_pred, k)

# A 5 años----
horizon <- 5

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk5y

# First calculate Aalen-Johansen estimate (as 'observed')
obj <- summary(survfit(Surv(time5y, status) ~ 1, 
                        data = vdata), 
                times = horizon)

aj <- list(
  "obs" = obj$pstate[, primary_event + 1],
  "se" = obj$std.err[, primary_event + 1]
)


# Calculate O/E
OE <- aj$obs / mean(vdata$pred)

# For the confidence interval we use method proposed in Debray et al. (2017) doi:10.1136/bmj.i6460
k <- 2
alpha <- 0.05
OE_summary <- cbind(
  "OE" = OE,
  "Lower .95" = exp(log(OE - qnorm(1 - alpha / 2) * aj$se / aj$obs)),
  "Upper .95" = exp(log(OE + qnorm(1 - alpha / 2) * aj$se / aj$obs))
)

OE_summary5a <- round(OE_summary, k)

avg_obs <- cbind(
  "avg_obs" = aj$obs * 100,
  "Lower .95" = 100 * (aj$obs - qnorm(1 - alpha / 2) * aj$se),
  "Upper .95" = 100 * (aj$obs + qnorm(1 - alpha / 2) * aj$se)
)

avg_pred <- cbind(
  "avg_pred" = mean(vdata$pred) * 100
)

avg_obs5a <- round(avg_obs, k)
avg_pred5a <- round(avg_pred, k)

# Seleccion del grupo: Stages 3b-4----
vdata <- dataB %>% 
  select(id, risk2y, risk5y, time5y, eventd5y, eventd, time) %>%  
  drop_na() %>%  
  mutate(status = factor(eventd5y, 
                         levels = c(0, 1, 2), 
                         labels = c("Alive w/o Kidney Failure", 
                                    "Kidney Failure", 
                                    "Death w/o Kidney Failure")))
primary_event <- 1

# A 2 años----
horizon <- 2

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk2y


# First calculate Aalen-Johansen estimate (as 'observed')
obj <- summary(survfit(Surv(time5y, status) ~ 1, 
                        data = vdata), 
                times = horizon)

aj <- list(
  "obs" = obj$pstate[, primary_event + 1],
  "se" = obj$std.err[, primary_event + 1]
)

# Calculate O/E
OE <- aj$obs / mean(vdata$pred)

# For the confidence interval we use method proposed in Debray et al. (2017) doi:10.1136/bmj.i6460
k <- 2
alpha <- 0.05
OE_summary <- cbind(
  "OE" = OE,
  "Lower .95" = exp(log(OE - qnorm(1 - alpha / 2) * aj$se / aj$obs)),
  "Upper .95" = exp(log(OE + qnorm(1 - alpha / 2) * aj$se / aj$obs))
)

OE_summary2b <- round(OE_summary, k)

avg_obs <- cbind(
  "avg_obs" = aj$obs * 100,
  "Lower .95" = 100 * (aj$obs - qnorm(1 - alpha / 2) * aj$se),
  "Upper .95" = 100 * (aj$obs + qnorm(1 - alpha / 2) * aj$se)
)

avg_pred <- cbind(
  "avg_pred" = mean(vdata$pred) * 100
)

avg_obs2b <- round(avg_obs, k)
avg_pred2b <- round(avg_pred, k)

# A 5 años----
horizon <- 5

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk5y


# First calculate Aalen-Johansen estimate (as 'observed')
obj <- summary(survfit(Surv(time5y, status) ~ 1, 
                        data = vdata), 
                times = horizon)

aj <- list(
  "obs" = obj$pstate[, primary_event + 1],
  "se" = obj$std.err[, primary_event + 1]
)


# Calculate O/E
OE <- aj$obs / mean(vdata$pred)

# For the confidence interval we use method proposed in Debray et al. (2017) doi:10.1136/bmj.i6460
k <- 2
alpha <- 0.05
OE_summary <- cbind(
  "OE" = OE,
  "Lower .95" = exp(log(OE - qnorm(1 - alpha / 2) * aj$se / aj$obs)),
  "Upper .95" = exp(log(OE + qnorm(1 - alpha / 2) * aj$se / aj$obs))
)

OE_summary5b <- round(OE_summary, k)

avg_obs <- cbind(
  "avg_obs" = aj$obs * 100,
  "Lower .95" = 100 * (aj$obs - qnorm(1 - alpha / 2) * aj$se),
  "Upper .95" = 100 * (aj$obs + qnorm(1 - alpha / 2) * aj$se)
)

avg_pred <- cbind(
  "avg_pred" = mean(vdata$pred) * 100
)

avg_obs5b <- round(avg_obs, k)
avg_pred5b <- round(avg_pred, k)

Weak calibration: Calibration intercept and Calibration slope

Mostrar código
# Seleccion del grupo: Stages 3-4----
vdata <- dataA %>% 
  select(id, risk2y, risk5y, time5y, eventd5y, eventd, time) %>%  
  drop_na() %>%  
  mutate(status = factor(eventd5y, 
                         levels = c(0, 1, 2), 
                         labels = c("Alive w/o Kidney Failure", 
                                    "Kidney Failure", 
                                    "Death w/o Kidney Failure")))
primary_event <- 1

# A 2 años----
horizon <- 2

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk2y

# First compute riskRegression::Score()
score_vdata <- Score(
  list("csh_validation" = vdata$pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  #  metrics = c("auc", "brier"),
  summary = c("ipa"),
  cause = primary_event,
  plots = "calibration"
)

# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)
pseudos <- data.frame(score_vdata$Calibration$plotframe)

# Note:
# - 'pseudos' is the data.frame with ACTUAL pseudo-observations, not the smoothed ones
# - Column ID is not the id in vdata; it is just a number assigned to each row of
# the original validation data sorted by time and event indicator
pseudos$cll_pred <- log(-log(1 - pseudos$risk)) # add the cloglog risk ests

# Fit model for calibration intercept
fit_cal_int <- geese(
  pseudovalue ~ offset(cll_pred),
  data = pseudos,
  id = ID,
  scale.fix = TRUE,
  family = gaussian,
  mean.link = "cloglog",
  corstr = "independence",
  jack = TRUE
)

# Fit model for calibration slope
fit_cal_slope <- geese(
  pseudovalue ~ offset(cll_pred) + cll_pred,
  data = pseudos,
  id = ID,
  scale.fix = TRUE,
  family = gaussian,
  mean.link = "cloglog",
  corstr = "independence",
  jack = TRUE
)

# Perform joint test on intercept and slope
betas <- fit_cal_slope$beta
vcov_mat <- fit_cal_slope$vbeta
wald <- drop(betas %*% solve(vcov_mat) %*% betas)
# pchisq(wald, df = 2, lower.tail = FALSE)

k <- 2
res_cal <- rbind(

  # Value, confidence interval and test for calibration intercept
  "Intercept, Stages 3-4, 2 year" = with(
    summary(fit_cal_int)$mean,
    c(
      "estimate" = estimate,
      `2.5 %` = estimate - qnorm(0.975) * san.se,
      `97.5 %` = estimate + qnorm(0.975) * san.se
    )
  ),

  # Value, confidence interval and test for calibration slope
  "Slope, Stages 3-4, 2 year" = with(
    summary(fit_cal_slope)$mean["cll_pred", ],
    c(
      "estimate" = 1 + estimate,
      `2.5 %` = 1 + (estimate - qnorm(0.975) * san.se),
      `97.5 %` = 1 + (estimate + qnorm(0.975) * san.se)
    )
  )
)

res_cal2a <- round(res_cal, k)

# A 5 años----
horizon <- 5

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk5y

# First compute riskRegression::Score()
score_vdata <- Score(
  list("csh_validation" = vdata$pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  #  metrics = c("auc", "brier"),
  summary = c("ipa"),
  cause = primary_event,
  plots = "calibration"
)

# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)
pseudos <- data.frame(score_vdata$Calibration$plotframe)

# Note:
# - 'pseudos' is the data.frame with ACTUAL pseudo-observations, not the smoothed ones
# - Column ID is not the id in vdata; it is just a number assigned to each row of
# the original validation data sorted by time and event indicator
pseudos$cll_pred <- log(-log(1 - pseudos$risk)) # add the cloglog risk ests

# Fit model for calibration intercept
fit_cal_int <- geese(
  pseudovalue ~ offset(cll_pred),
  data = pseudos,
  id = ID,
  scale.fix = TRUE,
  family = gaussian,
  mean.link = "cloglog",
  corstr = "independence",
  jack = TRUE
)

# Fit model for calibration slope
fit_cal_slope <- geese(
  pseudovalue ~ offset(cll_pred) + cll_pred,
  data = pseudos,
  id = ID,
  scale.fix = TRUE,
  family = gaussian,
  mean.link = "cloglog",
  corstr = "independence",
  jack = TRUE
)

# Perform joint test on intercept and slope
betas <- fit_cal_slope$beta
vcov_mat <- fit_cal_slope$vbeta
wald <- drop(betas %*% solve(vcov_mat) %*% betas)
# pchisq(wald, df = 2, lower.tail = FALSE)

k <- 2
res_cal <- rbind(

  # Value, confidence interval and test for calibration intercept
  "Intercept, Stages 3-4, 5 year" = with(
    summary(fit_cal_int)$mean,
    c(
      "estimate" = estimate,
      `2.5 %` = estimate - qnorm(0.975) * san.se,
      `97.5 %` = estimate + qnorm(0.975) * san.se
    )
  ),

  # Value, confidence interval and test for calibration slope
  "Slope, Stages 3-4, 5 year" = with(
    summary(fit_cal_slope)$mean["cll_pred", ],
    c(
      "estimate" = 1 + estimate,
      `2.5 %` = 1 + (estimate - qnorm(0.975) * san.se),
      `97.5 %` = 1 + (estimate + qnorm(0.975) * san.se)
    )
  )
)

res_cal5a <- round(res_cal, k)

# Seleccion del grupo: Stages 3b-4----

vdata <- dataB %>% 
  select(id, risk2y, risk5y, time5y, eventd5y, time, eventd) %>%  
  drop_na() %>%  
  mutate(status = factor(eventd5y, 
                         levels = c(0, 1, 2), 
                         labels = c("Alive w/o Kidney Failure", 
                                    "Kidney Failure", 
                                    "Death w/o Kidney Failure")))
primary_event <- 1

# A 2 años----
horizon <- 2

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk2y

# First compute riskRegression::Score()
score_vdata <- Score(
  list("csh_validation" = vdata$pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  #  metrics = c("auc", "brier"),
  summary = c("ipa"),
  cause = primary_event,
  plots = "calibration"
)

# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)
pseudos <- data.frame(score_vdata$Calibration$plotframe)

# Note:
# - 'pseudos' is the data.frame with ACTUAL pseudo-observations, not the smoothed ones
# - Column ID is not the id in vdata; it is just a number assigned to each row of
# the original validation data sorted by time and event indicator
pseudos$cll_pred <- log(-log(1 - pseudos$risk)) # add the cloglog risk ests

# Fit model for calibration intercept
fit_cal_int <- geese(
  pseudovalue ~ offset(cll_pred),
  data = pseudos,
  id = ID,
  scale.fix = TRUE,
  family = gaussian,
  mean.link = "cloglog",
  corstr = "independence",
  jack = TRUE
)

# Fit model for calibration slope
fit_cal_slope <- geese(
  pseudovalue ~ offset(cll_pred) + cll_pred,
  data = pseudos,
  id = ID,
  scale.fix = TRUE,
  family = gaussian,
  mean.link = "cloglog",
  corstr = "independence",
  jack = TRUE
)

# Perform joint test on intercept and slope
betas <- fit_cal_slope$beta
vcov_mat <- fit_cal_slope$vbeta
wald <- drop(betas %*% solve(vcov_mat) %*% betas)
# pchisq(wald, df = 2, lower.tail = FALSE)

k <- 2
res_cal <- rbind(

  # Value, confidence interval and test for calibration intercept
  "Intercept, Stages 3b-4, 2 year" = with(
    summary(fit_cal_int)$mean,
    c(
      "estimate" = estimate,
      `2.5 %` = estimate - qnorm(0.975) * san.se,
      `97.5 %` = estimate + qnorm(0.975) * san.se
    )
  ),

  # Value, confidence interval and test for calibration slope
  "Slope, Stages 3b-4, 2 year" = with(
    summary(fit_cal_slope)$mean["cll_pred", ],
    c(
      "estimate" = 1 + estimate,
      `2.5 %` = 1 + (estimate - qnorm(0.975) * san.se),
      `97.5 %` = 1 + (estimate + qnorm(0.975) * san.se)
    )
  )
)

res_cal2b <- round(res_cal, k)

# A 5 años----
horizon <- 5
vdata$pred <- vdata$risk5y

# First compute riskRegression::Score()
score_vdata <- Score(
  list("csh_validation" = vdata$pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  #  metrics = c("auc", "brier"),
  summary = c("ipa"),
  cause = primary_event,
  plots = "calibration"
)

# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)
pseudos <- data.frame(score_vdata$Calibration$plotframe)

# Note:
# - 'pseudos' is the data.frame with ACTUAL pseudo-observations, not the smoothed ones
# - Column ID is not the id in vdata; it is just a number assigned to each row of
# the original validation data sorted by time and event indicator
pseudos$cll_pred <- log(-log(1 - pseudos$risk)) # add the cloglog risk ests

# Fit model for calibration intercept
fit_cal_int <- geese(
  pseudovalue ~ offset(cll_pred),
  data = pseudos,
  id = ID,
  scale.fix = TRUE,
  family = gaussian,
  mean.link = "cloglog",
  corstr = "independence",
  jack = TRUE
)

# Fit model for calibration slope
fit_cal_slope <- geese(
  pseudovalue ~ offset(cll_pred) + cll_pred,
  data = pseudos,
  id = ID,
  scale.fix = TRUE,
  family = gaussian,
  mean.link = "cloglog",
  corstr = "independence",
  jack = TRUE
)

# Perform joint test on intercept and slope
betas <- fit_cal_slope$beta
vcov_mat <- fit_cal_slope$vbeta
wald <- drop(betas %*% solve(vcov_mat) %*% betas)
# pchisq(wald, df = 2, lower.tail = FALSE)

k <- 2
res_cal <- rbind(

  # Value, confidence interval and test for calibration intercept
  "Intercept, Stages 3b-4, 5 year" = with(
    summary(fit_cal_int)$mean,
    c(
      "estimate" = estimate,
      `2.5 %` = estimate - qnorm(0.975) * san.se,
      `97.5 %` = estimate + qnorm(0.975) * san.se
    )
  ),

  # Value, confidence interval and test for calibration slope
  "Slope, Stages 3b-4, 5 year" = with(
    summary(fit_cal_slope)$mean["cll_pred", ],
    c(
      "estimate" = 1 + estimate,
      `2.5 %` = 1 + (estimate - qnorm(0.975) * san.se),
      `97.5 %` = 1 + (estimate + qnorm(0.975) * san.se)
    )
  )
)

res_cal5b <- round(res_cal, k)
Calibration-in-the-large
  • OE ratio is a measure of calibration-in-the-large.
  • Calibration intercept is also a measure of calibration-in-the-large.

Moderate calibration: Calibration curves lowess based on pseudovalues

Mostrar código
# Seleccion del grupo: Stages 3-4----
vdata <- dataA %>% 
  select(id, risk2y, risk5y, time5y, eventd5y, time, eventd) %>%  
  drop_na()

primary_event <- 1

# A 2 años----
horizon <- 2

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk2y
pred <- as.matrix(vdata$pred)

# Calibration plot (pseudo-obs approach) ----------------------------------
# First compute riskRegression::Score()
score_vdata <- Score(
  list("csh_validation" = pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  #  metrics = c("auc", "brier"),
  summary = c("ipa"),
  cause = primary_event,
  plots = "calibration"
)

# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)
pseudos <- data.frame(score_vdata$Calibration$plotframe)
pseudos <- pseudos[order(pseudos$risk), ]

# Use linear loess (weighted local regression with polynomial degree = 1) smoothing
smooth_pseudos <- predict(
  stats::loess(pseudovalue ~ risk, data = pseudos, degree = 1, span = 0.33), 
  se = TRUE
)

pseudo_vals <- data.frame(
  obs = smooth_pseudos$fit, 
  risk = pseudos$risk, 
  se = smooth_pseudos$se, 
  df = smooth_pseudos$df
  ) %>% 
  mutate(
    ll = pmax(obs - qt(0.975, df) * se, 0), 
    ul = obs + qt(0.975, df) * se
  )

pseudo_vals %>% 
  ggplot(aes(x = risk, y = obs)) +
  geom_ribbon(aes(ymin = ll, ymax = ul), fill = "grey90") + 
  geom_abline(intercept = 0, slope = 1, colour = "red", linetype = 2) + 
  geom_line() + 
  scale_y_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) + 
  scale_x_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) + 
  theme_bw() + 
  xlab("Predicted risks") + 
  ylab("Observed outcome proportions") -> p1

# A 5 años----
horizon <- 5

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk5y
pred <- as.matrix(vdata$pred)

# Calibration plot (pseudo-obs approach) ----------------------------------
# First compute riskRegression::Score()
score_vdata <- Score(
  list("csh_validation" = pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  #  metrics = c("auc", "brier"),
  summary = c("ipa"),
  cause = primary_event,
  plots = "calibration"
)

# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)
pseudos <- data.frame(score_vdata$Calibration$plotframe)
pseudos <- pseudos[order(pseudos$risk), ]

# Use linear loess (weighted local regression with polynomial degree = 1) smoothing
smooth_pseudos <- predict(
  stats::loess(pseudovalue ~ risk, data = pseudos, degree = 1, span = 0.33), 
  se = TRUE
)

pseudo_vals <- data.frame(
  obs = smooth_pseudos$fit, 
  risk = pseudos$risk, 
  se = smooth_pseudos$se, 
  df = smooth_pseudos$df
  ) %>% 
  mutate(
    ll = pmax(obs - qt(0.975, df) * se, 0), 
    ul = obs + qt(0.975, df) * se
  )

pseudo_vals %>% 
  ggplot(aes(x = risk, y = obs)) +
  geom_ribbon(aes(ymin = ll, ymax = ul), fill = "grey90") + 
  geom_abline(intercept = 0, slope = 1, colour = "red", linetype = 2) + 
  geom_line() + 
  scale_y_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) + 
  scale_x_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) +   
  xlab("Predicted risks") + 
  ylab("Observed outcome proportions") + 
  theme_bw() -> p2

# Seleccion del grupo: Stages 3b-4----
vdata <- dataB %>% 
  select(id, risk2y, risk5y, time5y, eventd5y, time, eventd) %>%  
  drop_na()

primary_event <- 1

# A 2 años----
horizon <- 2

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk2y
pred <- as.matrix(vdata$pred)

# Calibration plot (pseudo-obs approach) ----------------------------------
# First compute riskRegression::Score()
score_vdata <- Score(
  list("csh_validation" = pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  #  metrics = c("auc", "brier"),
  summary = c("ipa"),
  cause = primary_event,
  plots = "calibration"
)

# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)
pseudos <- data.frame(score_vdata$Calibration$plotframe)
pseudos <- pseudos[order(pseudos$risk), ]

# Use linear loess (weighted local regression with polynomial degree = 1) smoothing
smooth_pseudos <- predict(
  stats::loess(pseudovalue ~ risk, data = pseudos, degree = 1, span = 0.33), 
  se = TRUE
)

pseudo_vals <- data.frame(
  obs = smooth_pseudos$fit, 
  risk = pseudos$risk, 
  se = smooth_pseudos$se, 
  df = smooth_pseudos$df
  ) %>% 
  mutate(
    ll = pmax(obs - qt(0.975, df) * se, 0), 
    ul = obs + qt(0.975, df) * se
  )

pseudo_vals %>% 
  ggplot(aes(x = risk, y = obs)) +
  geom_ribbon(aes(ymin = ll, ymax = ul), fill = "grey90") + 
  geom_abline(intercept = 0, slope = 1, colour = "red", linetype = 2) + 
  geom_line() + 
  scale_y_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) + 
  scale_x_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) + 
  theme_bw() + 
  xlab("Predicted risks") + 
  ylab("Observed outcome proportions") -> p3

# A 5 años----
horizon <- 5

# Add estimated risk and complementary log-log of it to dataset
vdata$pred <- vdata$risk5y
pred <- as.matrix(vdata$pred)

# Calibration plot (pseudo-obs approach) ----------------------------------
# First compute riskRegression::Score()
score_vdata <- Score(
  list("csh_validation" = pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  #  metrics = c("auc", "brier"),
  summary = c("ipa"),
  cause = primary_event,
  plots = "calibration"
)

# Use pseudo-observations calculated by Score() (can alternatively use pseudo::pseudoci)
pseudos <- data.frame(score_vdata$Calibration$plotframe)
pseudos <- pseudos[order(pseudos$risk), ]

# Use linear loess (weighted local regression with polynomial degree = 1) smoothing
smooth_pseudos <- predict(
  stats::loess(pseudovalue ~ risk, data = pseudos, degree = 1, span = 0.33), 
  se = TRUE
)

pseudo_vals <- data.frame(
  obs = smooth_pseudos$fit, 
  risk = pseudos$risk, 
  se = smooth_pseudos$se, 
  df = smooth_pseudos$df
  ) %>% 
  mutate(
    ll = pmax(obs - qt(0.975, df) * se, 0), 
    ul = obs + qt(0.975, df) * se
  )

pseudo_vals %>% 
  ggplot(aes(x = risk, y = obs)) +
  geom_ribbon(aes(ymin = ll, ymax = ul), fill = "grey90") + 
  geom_abline(intercept = 0, slope = 1, colour = "red", linetype = 2) + 
  geom_line() + 
  scale_y_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) + 
  scale_x_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) +   
  xlab("Predicted risks") + 
  ylab("Observed outcome proportions") + 
  theme_bw() -> p4

p1a <- p1 + 
  labs(title = "CKD Stages 3a-3b-4\n(2 year KFRE)") + 
  theme(plot.title = element_text(hjust = 0.5))
p2a <- p2 + labs(title = "CKD Stages 3a-3b-4\n(5 year KFRE)") + 
  theme(plot.title = element_text(hjust = 0.5))
p3a <- p3 + labs(title = "CKD Stages 3b-4\n(2 year KFRE)") + 
  theme(plot.title = element_text(hjust = 0.5))
p4a <- p4 + labs(title = "CKD Stages 3b-4\n(5 year KFRE)") + 
  theme(plot.title = element_text(hjust = 0.5))

(p1a | p2a) / (p3a | p4a) + plot_annotation(tag_levels = 'A') -> plot_calibration

ggsave(filename = "Plot_Calibration.png", 
       device = "png", 
       plot = plot_calibration, 
       path = here("Figures"), 
       scale = 2, 
       width = 2100, 
       height = 2100,
       units = "px", 
       dpi = 300)

5.2 Fig 3

Mostrar código
plot_calibration

Figure 2: Calibration plots for each group and prediction horizon. The predicted probability is shown on the x asis, and the observed kindey failure rate is given on the y axis

5.3 Discrimination

C-index and time-dependent C/D AUC

Mostrar código
B <- 1000

# Seleccion del grupo: Stages 3b-4----
vdata <- dataB %>% 
  select(id, risk2y, risk5y, time5y, eventd5y, time, eventd) %>%  
  drop_na() %>%  
  mutate(status = factor(eventd5y, 
                         levels = c(0, 1, 2), 
                         labels = c("Alive w/o Kidney Failure", 
                                    "Kidney Failure", 
                                    "Death w/o Kidney Failure")), 
         status_num = eventd5y)

primary_event <- 1

# A 2 años----
vdata$pred <- vdata$risk2y
horizon <- 2 
pred <- as.matrix(vdata$pred)

# Validation set
C_vdata <- pec::cindex(
  object = list("CauseSpecificCox" = pred),
  formula = Hist(time, eventd) ~ 1,
  cause = primary_event,
  eval.times = horizon,
  data = vdata
)$AppCindex$CauseSpecificCox

# Bootstraping C-index to calculate the bootstrap percentile confidence intervals
set.seed(777)

C_boot <- function(B, data) {
  cstat <- c()
  n <- nrow(data)
  for (i in 1:B) {
    ids <- sample(1:n, n, TRUE)
    data_boot <- data[ids, ]
    
    pec::cindex(
    object = list("CauseSpecificCox" = as.matrix(data_boot$pred)),
    formula = Hist(time, eventd) ~ 1,
    cause = primary_event,
    eval.times = horizon,
    data = data_boot
    )$AppCindex$CauseSpecificCox -> cstat[i]
  }
  return(cstat)
}
set.seed(777)
C_vdata_boot <- C_boot(B = B, data = vdata)

# Time-dependent AUC ---------

# Validation data
score_vdata <- Score(
  list("csh_validation" = pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  metrics = c("auc"),
  cause = primary_event,
  plots = "calibration"
)

alpha <- .05
k <- 3
res_discr_csh <- matrix(c(
  ## C-index
  # Validation CSH1
  C_vdata,
  quantile(C_vdata_boot, probs = alpha / 2),
  quantile(C_vdata_boot, probs = 1 - alpha / 2),

  ## Time-dependent AUC
  # Validation CSH1
  score_vdata$AUC$score$AUC,
  score_vdata$AUC$score$AUC - qnorm(1 - alpha / 2) * score_vdata$AUC$score$se,
  score_vdata$AUC$score$AUC + qnorm(1 - alpha / 2) * score_vdata$AUC$score$se
  ),
  nrow = 2, ncol = 3, byrow = T,
  dimnames = list(
  c("C-index, Stages 3b-4, 2 year", "Time dependent AUC, Stages 3b-4, 2 year"),
  c("Estimate", "Lower .95", "Upper .95")
  )
  )

res_discr_csh2b <- round(res_discr_csh, k)

# A 5 años----
vdata$pred <- vdata$risk5y
horizon <- 5 
pred <- as.matrix(vdata$pred)

# Validation set
C_vdata <- pec::cindex(
  object = list("CauseSpecificCox" = pred),
  formula = Hist(time, eventd) ~ 1,
  cause = primary_event,
  eval.times = horizon,
  data = vdata
)$AppCindex$CauseSpecificCox

# Bootstraping C-index to calculate the bootstrap percentile confidence intervals
set.seed(777)

C_boot <- function(B, data) {
  cstat <- c()
  n <- nrow(data)
  for (i in 1:B) {
    ids <- sample(1:n, n, TRUE)
    data_boot <- data[ids, ]
    
    pec::cindex(
    object = list("CauseSpecificCox" = as.matrix(data_boot$pred)),
    formula = Hist(time, eventd) ~ 1,
    cause = primary_event,
    eval.times = horizon,
    data = data_boot
    )$AppCindex$CauseSpecificCox -> cstat[i]
  }
  return(cstat)
}
set.seed(777)
C_vdata_boot <- C_boot(B = B, data = vdata)

# Time-dependent AUC ---------

# Validation data
score_vdata <- Score(
  list("csh_validation" = pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  metrics = c("auc"),
  cause = primary_event,
  plots = "calibration"
)

alpha <- .05
k <- 3
res_discr_csh <- matrix(c(
  ## C-index
  # Validation CSH1
  C_vdata,
  quantile(C_vdata_boot, probs = alpha / 2),
  quantile(C_vdata_boot, probs = 1 - alpha / 2),

  ## Time-dependent AUC
  # Validation CSH1
  score_vdata$AUC$score$AUC,
  score_vdata$AUC$score$AUC - qnorm(1 - alpha / 2) * score_vdata$AUC$score$se,
  score_vdata$AUC$score$AUC + qnorm(1 - alpha / 2) * score_vdata$AUC$score$se
  ),
  nrow = 2, ncol = 3, byrow = T,
  dimnames = list(
  c("C-index, Stages 3b-4, 5 year", "Time dependent AUC, Stages 3b-4, 5 year"),
  c("Estimate", "Lower .95", "Upper .95")
  )
  )

res_discr_csh5b <- round(res_discr_csh, k)
Mostrar código
# Seleccion del grupo: Stages 3a-3b-4----
vdata <- dataA %>% 
  select(id, risk2y, risk5y, time5y, eventd5y, time, eventd) %>%  
  drop_na() %>%  
  mutate(status = factor(eventd5y, 
                         levels = c(0, 1, 2), 
                         labels = c("Alive w/o Kidney Failure", 
                                    "Kidney Failure", 
                                    "Death w/o Kidney Failure")), 
         status_num = eventd5y)

primary_event <- 1

# A 2 años----
vdata$pred <- vdata$risk2y
horizon <- 2 
pred <- as.matrix(vdata$pred)

# Validation set
C_vdata <- pec::cindex(
  object = list("CauseSpecificCox" = pred),
  formula = Hist(time, eventd) ~ 1,
  cause = primary_event,
  eval.times = horizon,
  data = vdata
)$AppCindex$CauseSpecificCox

# Bootstraping C-index to calculate the bootstrap percentile confidence intervals
set.seed(777)

C_boot <- function(B, data) {
  cstat <- c()
  n <- nrow(data)
  for (i in 1:B) {
    ids <- sample(1:n, n, TRUE)
    data_boot <- data[ids, ]
    
    pec::cindex(
    object = list("CauseSpecificCox" = as.matrix(data_boot$pred)),
    formula = Hist(time, eventd) ~ 1,
    cause = primary_event,
    eval.times = horizon,
    data = data_boot
    )$AppCindex$CauseSpecificCox -> cstat[i]
  }
  return(cstat)
}

set.seed(777)
C_vdata_boot <- C_boot(B = B, data = vdata)

# Time-dependent AUC ---------

# Validation data
score_vdata <- Score(
  list("csh_validation" = pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  metrics = c("auc"),
  cause = primary_event,
  plots = "calibration"
)

alpha <- .05
k <- 3
res_discr_csh <- matrix(c(
  ## C-index
  # Validation CSH1
  C_vdata,
  quantile(C_vdata_boot, probs = alpha / 2),
  quantile(C_vdata_boot, probs = 1 - alpha / 2),

  ## Time-dependent AUC
  # Validation CSH1
  score_vdata$AUC$score$AUC,
  score_vdata$AUC$score$AUC - qnorm(1 - alpha / 2) * score_vdata$AUC$score$se,
  score_vdata$AUC$score$AUC + qnorm(1 - alpha / 2) * score_vdata$AUC$score$se
  ),
  nrow = 2, ncol = 3, byrow = T,
  dimnames = list(
  c("C-index, Stages 3-4, 2 year", "Time dependent AUC, Stages 3-4, 2 year"),
  c("Estimate", "Lower .95", "Upper .95")
  )
  )

res_discr_csh2a <- round(res_discr_csh, k)

# A 5 años----
vdata$pred <- vdata$risk5y
horizon <- 5 
pred <- as.matrix(vdata$pred)

# Validation set
C_vdata <- pec::cindex(
  object = list("CauseSpecificCox" = pred),
  formula = Hist(time, eventd) ~ 1,
  cause = primary_event,
  eval.times = horizon,
  data = vdata
)$AppCindex$CauseSpecificCox

# Bootstraping C-index to calculate the bootstrap percentile confidence intervals

set.seed(777)

C_boot <- function(B, data) {
  cstat <- c()
  n <- nrow(data)
  for (i in 1:B) {
    ids <- sample(1:n, n, TRUE)
    data_boot <- data[ids, ]
    
    pec::cindex(
    object = list("CauseSpecificCox" = as.matrix(data_boot$pred)),
    formula = Hist(time, eventd) ~ 1,
    cause = primary_event,
    eval.times = horizon,
    data = data_boot
    )$AppCindex$CauseSpecificCox -> cstat[i]
  }
  return(cstat)
}

set.seed(777)
C_vdata_boot <- C_boot(B = B, data = vdata)

# Time-dependent AUC ---------

# Validation data
score_vdata <- Score(
  list("csh_validation" = pred),
  formula = Hist(time, eventd) ~ 1,
  cens.model = "km",
  data = vdata,
  conf.int = TRUE,
  times = horizon,
  metrics = c("auc"),
  cause = primary_event,
  plots = "calibration"
)

alpha <- .05
k <- 3
res_discr_csh <- matrix(c(
  ## C-index
  # Validation CSH1
  C_vdata,
  quantile(C_vdata_boot, probs = alpha / 2),
  quantile(C_vdata_boot, probs = 1 - alpha / 2),

  ## Time-dependent AUC
  # Validation CSH1
  score_vdata$AUC$score$AUC,
  score_vdata$AUC$score$AUC - qnorm(1 - alpha / 2) * score_vdata$AUC$score$se,
  score_vdata$AUC$score$AUC + qnorm(1 - alpha / 2) * score_vdata$AUC$score$se
  ),
  nrow = 2, ncol = 3, byrow = T,
  dimnames = list(
  c("C-index, Stages 3-4, 5 year", "Time dependent AUC, Stages 3-4, 5 year"),
  c("Estimate", "Lower .95", "Upper .95")
  )
  )

res_discr_csh5a <- round(res_discr_csh, k)
Mostrar código
# Average predicted risk
avg_pred <- cbind("metrica" = "Average predicted risk", 
                 avg_pred2a, 
                 avg_pred5a, 
                 avg_pred2b, 
                 avg_pred5b)

colnames(avg_pred) <- c("metrica", "est2ya", "est5ya", "est2yb", "est5yb")

avg_pred %>% 
  as_tibble() %>% 
  mutate(
    est2ya = as.character(glue("{est2ya}%")), 
    est5ya = as.character(glue("{est5ya}%")),
    est2yb = as.character(glue("{est2yb}%")),
    est5yb = as.character(glue("{est5yb}%"))   
  ) -> avg_pred; avg_pred

# Average observed proportion
avg_obs <- cbind("metrica" = "Average observed proportion (95% CI)", 
                 avg_obs2a, 
                 avg_obs5a, 
                 avg_obs2b, 
                 avg_obs5b)

colnames(avg_obs) <- c("metrica", "OE2a", "ll2a", "ul2a", 
                          "OE5a", "ll5a", "ul5a", 
                          "OE2b", "ll2b", "ul2b", 
                          "OE5b", "ll5b", "ul5b")
avg_obs %>% 
  as_tibble() %>% 
  mutate(
    est2ya = as.character(glue("{OE2a}% ({ll2a}% to {ul2a}%)")), 
    est5ya = as.character(glue("{OE5a}% ({ll5a}% to {ul5a}%)")),
    est2yb = as.character(glue("{OE2b}% ({ll2b}% to {ul2b}%)")),
    est5yb = as.character(glue("{OE5b}% ({ll5b}% to {ul5b}%)"))   
  ) %>% 
  select(metrica, starts_with("est")) -> avg_obs; avg_obs

# OE summary
OE_summary <- cbind("metrica" = "O/E ratio (95% CI)", 
                    OE_summary2a, 
                    OE_summary5a, 
                    OE_summary2b, 
                    OE_summary5b); OE_summary

colnames(OE_summary) <- c("metrica", "OE2a", "ll2a", "ul2a", 
                         "OE5a", "ll5a", "ul5a", 
                         "OE2b", "ll2b", "ul2b", 
                         "OE5b", "ll5b", "ul5b")

OE_summary %>% 
  as_tibble() %>% 
  mutate(
    est2ya = as.character(glue("{OE2a} ({ll2a} to {ul2a})")), 
    est5ya = as.character(glue("{OE5a} ({ll5a} to {ul5a})")),
    est2yb = as.character(glue("{OE2b} ({ll2b} to {ul2b})")),
    est5yb = as.character(glue("{OE5b} ({ll5b} to {ul5b})"))   
    ) %>% 
  select(metrica, starts_with("est")) -> OE_summary; OE_summary

# Calibration slope e intercept
res_cal <- cbind("metrica" = c("Calibration intercept (95% CI)", " Calibration slope (95% CI)"), 
                 res_cal2a, 
                 res_cal5a, 
                 res_cal2b, 
                 res_cal5b); res_cal

colnames(res_cal) <- c("metrica", "OE2a", "ll2a", "ul2a", 
                       "OE5a", "ll5a", "ul5a", 
                       "OE2b", "ll2b", "ul2b", 
                       "OE5b", "ll5b", "ul5b")
res_cal <- 
  res_cal %>% 
  as_tibble() %>% 
  mutate(
    est2ya = as.character(glue("{OE2a} ({ll2a} to {ul2a})")), 
    est5ya = as.character(glue("{OE5a} ({ll5a} to {ul5a})")),
    est2yb = as.character(glue("{OE2b} ({ll2b} to {ul2b})")),
    est5yb = as.character(glue("{OE5b} ({ll5b} to {ul5b})"))   
  ) %>% 
  select(metrica, starts_with("est")) 

rbind(avg_pred, avg_obs, OE_summary, res_cal) -> table_performance; table_performance

# C-index y C/D AUC
res_discr_csh <- cbind("metrica" = c("C-index up to t-years (95% CI)", " C/D AUC, at t years (95% CI)"), 
                 res_discr_csh2a, 
                 res_discr_csh5a, 
                 res_discr_csh2b, 
                 res_discr_csh5b); res_discr_csh

colnames(res_discr_csh) <- c("metrica", "OE2a", "ll2a", "ul2a", 
                       "OE5a", "ll5a", "ul5a", 
                       "OE2b", "ll2b", "ul2b", 
                       "OE5b", "ll5b", "ul5b")
res_discr_csh <- 
  res_discr_csh %>% 
  as_tibble() %>% 
  mutate(
    est2ya = as.character(glue("{OE2a} ({ll2a} to {ul2a})")), 
    est5ya = as.character(glue("{OE5a} ({ll5a} to {ul5a})")),
    est2yb = as.character(glue("{OE2b} ({ll2b} to {ul2b})")),
    est5yb = as.character(glue("{OE5b} ({ll5b} to {ul5b})"))   
  ) %>% 
  select(metrica, starts_with("est")) 

rbind(avg_pred, avg_obs, OE_summary, res_cal, res_discr_csh) -> table_performance; table_performance
Mostrar código
table_performance %>% 
  mutate(grupo = c(rep("Calibration", 5), rep("Discrimination", 2))) %>% 
  relocate(grupo, .before = "metrica") %>%
  bind_rows() %>% 
  as_grouped_data(groups = "grupo") %>% 
  flextable::as_flextable(hide_grouplabel = TRUE)  %>%
  set_header_labels(
    metrica = "Validation aspect and performance measure", 
    est2ya = "t = 2 year", 
    est5ya = "t = 5 year", 
    est2yb = "t = 2 year", 
    est5yb = "t = 5 year" 
  ) %>% 
  add_header_row(
    values = c("Validation aspect and performance measure", "CKD Stages 3a-3b-4", "CKD Stages 3b-4"), 
    colwidths = c(1, 2, 2)
  ) %>% 
  merge_v(j = 1, part = "header") %>% 
  bold(i = c(1, 7)) %>% 
  autofit()  %>% 
  set_caption("Table 2. Performance measures of KFRE in the external dataset of patients with CKD Stages 3a-4 and 3b-4") %>% 
  theme_booktabs() %>%  
  bold(bold = TRUE, part = "header")  -> table_perf_final

table_perf_final %>% 
  flextable::save_as_docx(path = here("Tables/Table2.docx"))

5.4 Table 2

Mostrar código
table_perf_final

Validation aspect and performance measure

CKD Stages 3a-3b-4

CKD Stages 3b-4

t = 2 year

t = 5 year

t = 2 year

t = 5 year

Calibration

Average predicted risk

0.96%

3.18%

2.36%

7.66%

Average observed proportion (95% CI)

1.52% (1.24% to 1.79%)

3.37% (2.95% to 3.8%)

3.15% (2.5% to 3.79%)

6.86% (5.89% to 7.84%)

O/E ratio (95% CI)

1.57 (1.39 to 1.76)

1.06 (0.93 to 1.19)

1.33 (1.13 to 1.54)

0.9 (0.75 to 1.04)

Calibration intercept (95% CI)

0.18 (-0.1 to 0.45)

-0.26 (-0.45 to -0.07)

0.16 (-0.12 to 0.44)

-0.29 (-0.48 to -0.1)

Calibration slope (95% CI)

0.79 (0.61 to 0.96)

0.75 (0.65 to 0.86)

0.82 (0.6 to 1.03)

0.79 (0.66 to 0.92)

Discrimination

C-index up to t-years (95% CI)

0.853 (0.812 to 0.892)

0.845 (0.818 to 0.872)

0.848 (0.803 to 0.885)

0.827 (0.796 to 0.857)

C/D AUC, at t years (95% CI)

0.855 (0.816 to 0.895)

0.847 (0.819 to 0.875)

0.853 (0.811 to 0.895)

0.836 (0.803 to 0.869)

6 Supplementary tables

6.1 Table S1

Mostrar código
table_kfre <- data.frame(
  pred = c("2-years", "5-years"), 
  eq = c("$1-{0.9832}^{e^{(-0.2201\times(\frac{age}{10}-7.036)+0.2467\times(male-0.5642)-0.5567\times(\frac{eGFR}{5}-7.222)+0.4510\times(log{(ACR)}-5.137))}}$", 
         "$1-{0.9365}^{e^{(-0.2201\times(\frac{age}{10}-7.036)+0.2467\times(male-0.5642)-0.5567\times(\frac{eGFR}{5}-7.222)+0.4510\times(log{(ACR)}-5.137))}}$")
)
Mostrar código
knitr::kable(table_kfre, escape = TRUE, 
             col.names = c("Prediction horizons", 
                           "Original regional equation calibrated for predicted risk of kidney failure"), 
             caption = "**Table S1.** KFRE equations externally validated by the study")
**Table S1.** KFRE equations externally validated by the study
Prediction horizons Original regional equation calibrated for predicted risk of kidney failure
2-years $1-{0.9832}^{e^{(-0.2201 imes(rac{age}{10}-7.036)+0.2467 imes(male-0.5642)-0.5567 imes(rac{eGFR}{5}-7.222)+0.4510 imes(log{(ACR)}-5.137))}}$
5-years $1-{0.9365}^{e^{(-0.2201 imes(rac{age}{10}-7.036)+0.2467 imes(male-0.5642)-0.5567 imes(rac{eGFR}{5}-7.222)+0.4510 imes(log{(ACR)}-5.137))}}$

6.2 Table S2

Mostrar código
table_coding <- data.frame(
  Variable = c("age", "male", "eGFR_ckdepi", "acr"), 
  Coding = c("integer number that indicates the age in completed years", 
             "1 = male; 0 = female", 
             "estimated glomerular filtration rate obtained by CKD-EPI formula in $ml/min/1.73m^2$", 
             "albumin-to-creatinine ratio in mg/g")
)
Mostrar código
knitr::kable(table_coding, escape = TRUE, 
             caption = "*Table S2.* Coding of variables")
*Table S2.* Coding of variables
Variable Coding
age integer number that indicates the age in completed years
male 1 = male; 0 = female
eGFR_ckdepi estimated glomerular filtration rate obtained by CKD-EPI formula in $ml/min/1.73m^2$
acr albumin-to-creatinine ratio in mg/g

6.3 Table S3

Mostrar código
data %>%  
  filter(ckd_stage == "Stages 3-4") %>% 
  mutate(
    primary_event = case_when(
      eventd5ylab == "Kidney Failure" ~ "Kidney failure", 
      eventd5ylab %in% c("Alive w/o Kidney Failure", "Death w/o Kidney Failure") 
      ~ "No kidney failure"), 
    primary_event = factor(primary_event, 
                           levels = c("No kidney failure", 
                                      "Kidney failure"))) %>% 
  dplyr::select(sex, age, hta, dm, acr_cat, grf_cat, ckd_class, crea,
                eGFR_ckdepi, acr, urine_album, urine_crea,  death2y, 
                eventd2ylab, death5y, eventd5ylab, primary_event) %>% 
  gtsummary::tbl_summary(
  by = "primary_event", 
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({p25}, {p75})", 
    "{min}, {max}"
  ),
  digits = list(all_continuous() ~ c(1, 1, 1, 1), 
                all_categorical() ~ c(0, 1)), 
  missing_text = "Missing"
  ) %>% 
  bold_labels() -> tab2


data %>%  
  filter(ckd_stage == "Stages 3-4") %>% 
  mutate(
    primary_event = case_when(
      eventd2ylab == "Kidney Failure" ~ "Kidney failure", 
      eventd2ylab %in% c("Alive w/o Kidney Failure", "Death w/o Kidney Failure") 
      ~ "No kidney failure"), 
    primary_event = factor(primary_event, 
                           levels = c("No kidney failure", 
                                      "Kidney failure"))) %>% 
  dplyr::select(sex, age, hta, dm, acr_cat, grf_cat, ckd_class, crea,
                eGFR_ckdepi, acr, urine_album, urine_crea,  death2y, 
                eventd2ylab, death5y, eventd5ylab,   
                primary_event) %>% 
  gtsummary::tbl_summary(
  by = "primary_event", 
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({p25}, {p75})", 
    "{min}, {max}"
  ),
  digits = list(all_continuous() ~ c(1, 1, 1, 1), 
                all_categorical() ~ c(0, 1)), 
  missing_text = "Missing"
  ) %>% 
  bold_labels() -> tab1

tbl_merge(list(tab1, tab2), 
          tab_spanner = c("2-years", "5-years")) %>% 
  modify_caption("Table S3. Baseline characteristics of the study patients at CKD Stages 3a, 3b or 4 for all predictors, stratified by outcome at 2- and 5-years") -> tab


tab %>% 
  as_flex_table() %>% 
  flextable::save_as_docx(path = here("Tables/TableS3.docx"))

tab %>% 
  gtsummary::as_kable_extra() %>%
  kableExtra::kable_styling("striped")
Table S3. Baseline characteristics of the study patients at CKD Stages 3a, 3b or 4 for all predictors, stratified by outcome at 2- and 5-years
2-years
5-years
Characteristic No kidney failure, N = 7,405 Kidney failure, N = 114 No kidney failure, N = 7,280 Kidney failure, N = 239
Sex
Female 4,062 (54.9%) 45 (39.5%) 4,015 (55.2%) 92 (38.5%)
Male 3,343 (45.1%) 69 (60.5%) 3,265 (44.8%) 147 (61.5%)
Age (years)
Mean (SD) 74.1 (10.2) 66.6 (11.7) 74.2 (10.1) 66.5 (12.5)
Median (IQR) 75.0 (68.0, 81.0) 67.0 (59.2, 74.0) 75.0 (68.0, 82.0) 67.0 (59.0, 75.0)
Range 23.0, 97.0 36.0, 88.0 23.0, 97.0 26.0, 94.0
Hypertension 4,421 (59.7%) 65 (57.0%) 4,339 (59.6%) 147 (61.5%)
Diabetes Mellitus 1,796 (24.3%) 49 (43.0%) 1,743 (23.9%) 102 (42.7%)
Persistent albuminuria categories
A1 4,749 (64.1%) 23 (20.2%) 4,725 (64.9%) 47 (19.7%)
A2 1,984 (26.8%) 34 (29.8%) 1,941 (26.7%) 77 (32.2%)
A3 672 (9.1%) 57 (50.0%) 614 (8.4%) 115 (48.1%)
GFR categories
G3a 4,695 (63.4%) 26 (22.8%) 4,664 (64.1%) 57 (23.8%)
G3b 2,171 (29.3%) 36 (31.6%) 2,126 (29.2%) 81 (33.9%)
G4 539 (7.3%) 52 (45.6%) 490 (6.7%) 101 (42.3%)
CKD KDIGO classification
Moderately increased risk 3,266 (44.1%) 12 (10.5%) 3,256 (44.7%) 22 (9.2%)
High risk 2,448 (33.1%) 12 (10.5%) 2,429 (33.4%) 31 (13.0%)
Very high risk 1,691 (22.8%) 90 (78.9%) 1,595 (21.9%) 186 (77.8%)
Serum Creatinine (mg/dL)
Mean (SD) 1.4 (0.4) 2.0 (0.7) 1.3 (0.4) 2.0 (0.7)
Median (IQR) 1.3 (1.1, 1.5) 1.9 (1.5, 2.5) 1.3 (1.1, 1.5) 1.9 (1.5, 2.4)
Range 0.8, 3.8 1.0, 3.9 0.8, 3.8 0.9, 3.9
eGFR using CKD-EPI, ml/min/1.73m2
Mean (SD) 46.4 (9.6) 33.7 (12.5) 46.6 (9.4) 34.5 (12.4)
Median (IQR) 48.9 (40.8, 53.9) 31.7 (23.2, 42.8) 49.0 (41.1, 53.9) 33.0 (23.5, 44.2)
Range 15.0, 60.0 15.4, 59.8 15.1, 60.0 15.0, 59.8
Albumin-to-creatinine ratio, mg/g
Mean (SD) 235.5 (3,059.4) 1,101.7 (1,614.1) 229.3 (3,083.4) 836.1 (1,283.1)
Median (IQR) 14.2 (4.4, 62.7) 302.0 (52.5, 1,663.9) 14.0 (4.4, 59.4) 270.9 (51.2, 992.4)
Range 0.0, 144,870.6 2.5, 7,462.7 0.0, 144,870.6 0.2, 7,462.7
Urine Albumin (mg/ml)
Mean (SD) 7.6 (26.5) 52.1 (73.5) 7.1 (25.7) 44.0 (62.3)
Median (IQR) 0.9 (0.3, 3.7) 16.0 (4.3, 71.6) 0.9 (0.3, 3.5) 15.3 (3.4, 56.5)
Range 0.0, 658.0 0.2, 348.1 0.0, 658.0 0.0, 348.1
Urine Creatinine (mg/dl)
Mean (SD) 72.5 (47.7) 64.8 (33.7) 72.4 (47.0) 69.8 (61.0)
Median (IQR) 63.4 (41.4, 86.5) 59.6 (39.6, 85.0) 63.5 (41.3, 86.6) 58.7 (41.0, 85.0)
Range 0.1, 722.1 6.4, 218.6 0.1, 620.1 6.4, 722.1
Death at 2 years 563 (7.6%) 77 (67.5%) 563 (7.7%) 77 (32.2%)
Outcome at 2 years
Alive w/o Kidney Failure 6,842 (92.4%) 0 (0.0%) 6,717 (92.3%) 125 (52.3%)
Death w/o Kidney Failure 563 (7.6%) 0 (0.0%) 563 (7.7%) 0 (0.0%)
Kidney Failure 0 (0.0%) 114 (100.0%) 0 (0.0%) 114 (47.7%)
Death at 5 years 1,462 (19.7%) 77 (67.5%) 1,400 (19.2%) 139 (58.2%)
Outcome at 5 years
Alive w/o Kidney Failure 5,880 (79.4%) 0 (0.0%) 5,880 (80.8%) 0 (0.0%)
Death w/o Kidney Failure 1,400 (18.9%) 0 (0.0%) 1,400 (19.2%) 0 (0.0%)
Kidney Failure 125 (1.7%) 114 (100.0%) 0 (0.0%) 239 (100.0%)
1 n (%)

6.4 Table S4

Mostrar código
data %>%  
  filter(ckd_stage2 == "Stages 3b-4") %>% 
  mutate(
    primary_event = case_when(
      eventd5ylab == "Kidney Failure" ~ "Kidney failure", 
      eventd5ylab %in% c("Alive w/o Kidney Failure", "Death w/o Kidney Failure") 
      ~ "No kidney failure"), 
    primary_event = factor(primary_event, 
                           levels = c("No kidney failure", 
                                      "Kidney failure"))) %>% 
  dplyr::select(sex, age, hta, dm, acr_cat, grf_cat, ckd_class, crea,
                eGFR_ckdepi, acr, urine_album, urine_crea,  death2y, 
                eventd2ylab, death5y, eventd5ylab,   
                primary_event) %>% 
  gtsummary::tbl_summary(
  by = "primary_event", 
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({p25}, {p75})", 
    "{min}, {max}"
  ),
  digits = list(all_continuous() ~ c(1, 1, 1, 1), 
                all_categorical() ~ c(0, 1)), 
  missing_text = "Missing"
  ) %>% 
  bold_labels() -> tab2


data %>%  
  filter(ckd_stage2 == "Stages 3b-4") %>% 
  mutate(
    primary_event = case_when(
      eventd2ylab == "Kidney Failure" ~ "Kidney failure", 
      eventd2ylab %in% c("Alive w/o Kidney Failure", "Death w/o Kidney Failure") 
      ~ "No kidney failure"), 
    primary_event = factor(primary_event, 
                           levels = c("No kidney failure", 
                                      "Kidney failure"))) %>% 
  dplyr::select(sex, age, hta, dm, acr_cat, grf_cat, ckd_class, crea,
                eGFR_ckdepi, acr, urine_album, urine_crea,  death2y, 
                eventd2ylab, death5y, eventd5ylab,  
                primary_event) %>% 
  gtsummary::tbl_summary(
  by = "primary_event", 
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({p25}, {p75})", 
    "{min}, {max}"
  ),
  digits = list(all_continuous() ~ c(1, 1, 1, 1), 
                all_categorical() ~ c(0, 1)), 
  missing_text = "Missing"
  ) %>% 
  bold_labels() -> tab1

tbl_merge(list(tab1, tab2), 
          tab_spanner = c("2-years", "5-years")) %>% 
  modify_caption("Table S4. Baseline characteristics of the study patients at CKD Stages 3b or 4 for all predictors, stratified by outcome at 2- and 5-years") -> tab

tab %>% 
  as_flex_table() %>% 
  flextable::save_as_docx(path = here("Tables/TableS4.docx"))

tab %>% 
  gtsummary::as_kable_extra() %>%
  kableExtra::kable_styling("striped")
Table S4. Baseline characteristics of the study patients at CKD Stages 3b or 4 for all predictors, stratified by outcome at 2- and 5-years
2-years
5-years
Characteristic No kidney failure, N = 2,710 Kidney failure, N = 88 No kidney failure, N = 2,616 Kidney failure, N = 182
Sex
Female 1,363 (50.3%) 35 (39.8%) 1,329 (50.8%) 69 (37.9%)
Male 1,347 (49.7%) 53 (60.2%) 1,287 (49.2%) 113 (62.1%)
Age (years)
Mean (SD) 75.9 (10.4) 66.1 (12.1) 76.2 (10.2) 66.9 (12.4)
Median (IQR) 77.0 (70.0, 83.0) 67.0 (59.0, 74.0) 77.0 (70.0, 83.0) 67.0 (59.2, 75.0)
Range 23.0, 97.0 36.0, 88.0 23.0, 97.0 26.0, 94.0
Hypertension 1,581 (58.3%) 55 (62.5%) 1,517 (58.0%) 119 (65.4%)
Diabetes Mellitus 635 (23.4%) 39 (44.3%) 596 (22.8%) 78 (42.9%)
Persistent albuminuria categories
A1 1,483 (54.7%) 11 (12.5%) 1,469 (56.2%) 25 (13.7%)
A2 830 (30.6%) 30 (34.1%) 798 (30.5%) 62 (34.1%)
A3 397 (14.6%) 47 (53.4%) 349 (13.3%) 95 (52.2%)
GFR categories
G3a 0 (0.0%) 0 (0.0%) 0 (0.0%) 0 (0.0%)
G3b 2,171 (80.1%) 36 (40.9%) 2,126 (81.3%) 81 (44.5%)
G4 539 (19.9%) 52 (59.1%) 490 (18.7%) 101 (55.5%)
CKD KDIGO classification
Moderately increased risk 0 (0.0%) 0 (0.0%) 0 (0.0%) 0 (0.0%)
High risk 1,294 (47.7%) 8 (9.1%) 1,286 (49.2%) 16 (8.8%)
Very high risk 1,416 (52.3%) 80 (90.9%) 1,330 (50.8%) 166 (91.2%)
Serum Creatinine (mg/dL)
Mean (SD) 1.7 (0.4) 2.3 (0.6) 1.7 (0.4) 2.2 (0.6)
Median (IQR) 1.6 (1.4, 1.9) 2.1 (1.8, 2.5) 1.6 (1.4, 1.8) 2.1 (1.7, 2.5)
Range 1.1, 3.8 1.2, 3.9 1.1, 3.8 1.2, 3.9
eGFR using CKD-EPI, ml/min/1.73m2
Mean (SD) 35.9 (7.1) 28.2 (8.1) 36.1 (7.0) 29.0 (8.3)
Median (IQR) 37.6 (31.8, 41.8) 26.9 (21.6, 34.9) 37.8 (32.1, 41.9) 28.3 (22.4, 35.4)
Range 15.0, 45.0 15.4, 43.8 15.1, 45.0 15.0, 44.9
Albumin-to-creatinine ratio, mg/g
Mean (SD) 306.8 (3,082.4) 1,172.8 (1,631.1) 296.0 (3,133.9) 881.4 (1,272.2)
Median (IQR) 24.6 (6.3, 145.4) 367.7 (149.1, 1,811.8) 23.3 (6.1, 131.2) 334.6 (143.8, 1,076.5)
Range 0.0, 144,870.6 2.5, 7,462.7 0.0, 144,870.6 1.4, 7,462.7
Urine Albumin (mg/ml)
Mean (SD) 11.7 (33.3) 56.1 (77.1) 10.7 (32.2) 47.4 (64.4)
Median (IQR) 1.5 (0.4, 8.7) 16.3 (8.9, 66.6) 1.4 (0.4, 7.7) 16.5 (8.6, 60.3)
Range 0.0, 658.0 0.2, 348.1 0.0, 658.0 0.1, 348.1
Urine Creatinine (mg/dl)
Mean (SD) 71.7 (47.4) 62.1 (33.5) 71.7 (45.8) 68.3 (61.8)
Median (IQR) 65.1 (43.7, 85.0) 55.3 (38.1, 85.0) 65.4 (43.9, 85.0) 58.0 (40.1, 85.0)
Range 0.7, 722.1 6.4, 218.6 0.7, 620.1 6.4, 722.1
Death at 2 years 300 (11.1%) 58 (65.9%) 300 (11.5%) 58 (31.9%)
Outcome at 2 years
Alive w/o Kidney Failure 2,410 (88.9%) 0 (0.0%) 2,316 (88.5%) 94 (51.6%)
Death w/o Kidney Failure 300 (11.1%) 0 (0.0%) 300 (11.5%) 0 (0.0%)
Kidney Failure 0 (0.0%) 88 (100.0%) 0 (0.0%) 88 (48.4%)
Death at 5 years 726 (26.8%) 58 (65.9%) 683 (26.1%) 101 (55.5%)
Outcome at 5 years
Alive w/o Kidney Failure 1,933 (71.3%) 0 (0.0%) 1,933 (73.9%) 0 (0.0%)
Death w/o Kidney Failure 683 (25.2%) 0 (0.0%) 683 (26.1%) 0 (0.0%)
Kidney Failure 94 (3.5%) 88 (100.0%) 0 (0.0%) 182 (100.0%)
1 n (%)

6.5 Table S5

Mostrar código
data %>%  
  dplyr::select(sex, age, hta, dm, acr_cat, grf_cat, ckd_class, crea,
                eGFR_ckdepi, acr, urine_album, urine_crea,  death2y, 
                eventd2ylab, death5y, eventd5ylab) %>% 
  gtsummary::tbl_summary(
  by = grf_cat, 
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({p25} - {p75})", 
    "{min} - {max}"
  ),
  digits = list(all_continuous() ~ c(1, 1, 1, 1), 
                all_categorical() ~ c(0, 1)), 
  missing_text = "Missing"
  ) %>% 
  # add_p() %>% 
  bold_labels() %>% 
  modify_caption("Table S5. Baseline characteristics and frequency of outcomes of the study population according to CKD stages") -> tab

tab %>% 
  as_flex_table() %>% 
  flextable::save_as_docx(path = here("Tables/TableS5.docx"))

tab %>% 
  gtsummary::as_kable_extra() %>%
  kableExtra::kable_styling("striped")
Table S5. Baseline characteristics and frequency of outcomes of the study population according to CKD stages
Characteristic G3a, N = 4,721 G3b, N = 2,207 G4, N = 591
Sex
Female 2,709 (57.4%) 1,102 (49.9%) 296 (50.1%)
Male 2,012 (42.6%) 1,105 (50.1%) 295 (49.9%)
Age (years)
Mean (SD) 73.1 (9.9) 76.0 (10.2) 73.9 (11.8)
Median (IQR) 74.0 (67.0 - 80.0) 77.0 (70.0 - 83.0) 75.0 (67.0 - 83.0)
Range 23.0 - 97.0 26.0 - 97.0 23.0 - 95.0
Hypertension 2,850 (60.4%) 1,295 (58.7%) 341 (57.7%)
Diabetes Mellitus 1,171 (24.8%) 521 (23.6%) 153 (25.9%)
Persistent albuminuria categories
A1 3,278 (69.4%) 1,302 (59.0%) 192 (32.5%)
A2 1,158 (24.5%) 632 (28.6%) 228 (38.6%)
A3 285 (6.0%) 273 (12.4%) 171 (28.9%)
CKD KDIGO classification
Moderately increased risk 3,278 (69.4%) 0 (0.0%) 0 (0.0%)
High risk 1,158 (24.5%) 1,302 (59.0%) 0 (0.0%)
Very high risk 285 (6.0%) 905 (41.0%) 591 (100.0%)
Serum Creatinine (mg/dL)
Mean (SD) 1.2 (0.2) 1.5 (0.2) 2.3 (0.5)
Median (IQR) 1.1 (1.0 - 1.3) 1.5 (1.3 - 1.7) 2.3 (1.9 - 2.6)
Range 0.8 - 1.9 1.1 - 2.4 1.5 - 3.9
eGFR using CKD-EPI, ml/min/1.73m2
Mean (SD) 52.5 (3.9) 38.8 (4.2) 24.2 (4.1)
Median (IQR) 52.6 (49.4 - 55.6) 39.3 (35.3 - 42.4) 24.7 (21.2 - 27.7)
Range 45.0 - 60.0 30.0 - 45.0 15.0 - 30.0
Albumin-to-creatinine ratio, mg/g
Mean (SD) 198.0 (3,039.8) 302.6 (3,386.3) 451.7 (1,109.7)
Median (IQR) 10.8 (3.9 - 41.5) 20.3 (5.3 - 100.1) 127.2 (17.7 - 405.2)
Range 0.0 - 137,672.1 0.0 - 144,870.6 0.1 - 18,259.2
Urine Albumin (mg/ml)
Mean (SD) 5.5 (21.8) 10.0 (31.4) 25.0 (48.8)
Median (IQR) 0.7 (0.2 - 2.4) 1.3 (0.3 - 5.8) 7.4 (1.3 - 20.9)
Range 0.0 - 548.0 0.0 - 658.0 0.0 - 383.0
Urine Creatinine (mg/dl)
Mean (SD) 72.9 (47.9) 71.3 (46.0) 72.0 (50.7)
Median (IQR) 62.4 (40.4 - 88.1) 64.8 (43.0 - 85.0) 65.3 (44.6 - 85.0)
Range 0.1 - 438.3 0.8 - 620.1 0.7 - 722.1
Death at 2 years 282 (6.0%) 230 (10.4%) 128 (21.7%)
Outcome at 2 years
Alive w/o Kidney Failure 4,432 (93.9%) 1,965 (89.0%) 445 (75.3%)
Death w/o Kidney Failure 263 (5.6%) 206 (9.3%) 94 (15.9%)
Kidney Failure 26 (0.6%) 36 (1.6%) 52 (8.8%)
Death at 5 years 755 (16.0%) 555 (25.1%) 229 (38.7%)
Outcome at 5 years
Alive w/o Kidney Failure 3,947 (83.6%) 1,622 (73.5%) 311 (52.6%)
Death w/o Kidney Failure 717 (15.2%) 504 (22.8%) 179 (30.3%)
Kidney Failure 57 (1.2%) 81 (3.7%) 101 (17.1%)
1 n (%)

6.6 Table S6

Mostrar código
data %>%  
  filter(ckd_stage == "Stages 3-4") %>% 
  dplyr::select(cas) %>% 
  gtsummary::tbl_summary(
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({min}, {max})"
  ),
  digits = list(all_categorical() ~ c(0, 1)), 
  percent = "row"
  ) %>%
  bold_labels() -> tab1a

data %>%  
  filter(ckd_stage2 == "Stages 3b-4") %>% 
  dplyr::select(cas) %>% 
  gtsummary::tbl_summary(
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({min}, {max})"
  ),
  digits = list(all_categorical() ~ c(0, 1)), 
  percent = "row"
  ) %>%
  bold_labels() -> tab1b

data %>%  
  filter(ckd_stage == "Stages 3-4") %>% 
  dplyr::select(cas, eventd5ylab) %>%  
  gtsummary::tbl_summary(
  by = "eventd5ylab",   
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({min}, {max})"
  ),
  digits = list(all_categorical() ~ c(0, 1)), 
  percent = "row"
  ) %>%  
  bold_labels() -> tab2a

data %>%  
  filter(ckd_stage2 == "Stages 3b-4") %>% 
  dplyr::select(cas, eventd5ylab) %>%  
  gtsummary::tbl_summary(
  by = "eventd5ylab",   
  type = all_continuous() ~ "continuous2",
  statistic = all_continuous() ~ c(
    "{mean} ({sd})",
    "{median} ({min}, {max})"
  ),
  digits = list(all_categorical() ~ c(0, 1)), 
  percent = "row"
  ) %>%  
  bold_labels() -> tab2b

tbl_merge(list(tab1a, tab2a), 
          tab_spanner = c("Overall", "Outcome")) %>% 
  modify_caption("Table S6. Distribution of CKD 3a-4 patients included in the analysis within 17 health facilities of the EsSalud Rebagliati Network") -> tabA

tabA %>% 
  as_flex_table() %>% 
  theme_booktabs() %>%  
  bold(bold = TRUE, part = "header") %>% 
  flextable::save_as_docx(path = here("Tables/TableS6.docx"))
Mostrar código
tabA %>% 
  gtsummary::as_kable_extra() %>% 
  add_header_above(c(" ", "CKD Stages 3a-3b-4" = 4)) %>% 
  kableExtra::kable_styling("striped")
Table S6. Distribution of CKD 3a-4 patients included in the analysis within 17 health facilities of the EsSalud Rebagliati Network
CKD Stages 3a-3b-4
Overall
Outcome
Characteristic N = 7,519 Alive w/o Kidney Failure, N = 5,880 Death w/o Kidney Failure, N = 1,400 Kidney Failure, N = 239
Healthcare center
C. M. MALA 74 (100.0%) 48 (64.9%) 18 (24.3%) 8 (10.8%)
CAP II LURIN 58 (100.0%) 52 (89.7%) 6 (10.3%) 0 (0.0%)
CAP III LOS PROCEDES DE SJM 403 (100.0%) 319 (79.2%) 73 (18.1%) 11 (2.7%)
CAP III SAN ISIDRO 287 (100.0%) 204 (71.1%) 71 (24.7%) 12 (4.2%)
CAP III SURQUILLO 308 (100.0%) 242 (78.6%) 53 (17.2%) 13 (4.2%)
HOSP. II ANGAMOS 1,471 (100.0%) 1,088 (74.0%) 298 (20.3%) 85 (5.8%)
HOSPITAL EDGARDO REBAGLIATI MARTINS 34 (100.0%) 20 (58.8%) 11 (32.4%) 3 (8.8%)
HOSPITAL I CARLOS ALCANTARA B. 1,032 (100.0%) 858 (83.1%) 160 (15.5%) 14 (1.4%)
HOSPITAL I ULDARICO ROCCA F. 207 (100.0%) 158 (76.3%) 42 (20.3%) 7 (3.4%)
HOSPITAL II CANETE 258 (100.0%) 179 (69.4%) 73 (28.3%) 6 (2.3%)
POLICLINICO CHEQUEOS LARCO 77 (100.0%) 71 (92.2%) 5 (6.5%) 1 (1.3%)
POLICLINICO CHINCHA 1,289 (100.0%) 1,058 (82.1%) 202 (15.7%) 29 (2.2%)
POLICLINICO JUAN JOSE RODRIGUEZ LAZO 679 (100.0%) 497 (73.2%) 166 (24.4%) 16 (2.4%)
POLICLINICO PABLO BERMUDEZ 752 (100.0%) 586 (77.9%) 143 (19.0%) 23 (3.1%)
POLICLINICO PROCERES 454 (100.0%) 381 (83.9%) 64 (14.1%) 9 (2.0%)
POLICLINICO SANTA CRUZ 102 (100.0%) 90 (88.2%) 10 (9.8%) 2 (2.0%)
POLICLINICO VILLA MARIA 34 (100.0%) 29 (85.3%) 5 (14.7%) 0 (0.0%)
1 n (%)

6.7 Table S7

Mostrar código
tbl_merge(list(tab1b, tab2b), 
          tab_spanner = c("Overall", "Outcome")) %>% 
  modify_caption("Table S7. Distribution of CKD 3b-4 patients included in the analysis within 17 health facilities of the EsSalud Rebagliati Network") -> tabB

tabB %>% 
  as_flex_table() %>% 
  theme_booktabs() %>%  
  bold(bold = TRUE, part = "header") %>% 
  flextable::save_as_docx(path = here("Tables/TableS7.docx"))
Mostrar código
tabB %>% 
  gtsummary::as_kable_extra() %>% 
  add_header_above(c(" ", "CKD Stages 3b-4" = 4)) %>% 
  kableExtra::kable_styling("striped")
Table S7. Distribution of CKD 3b-4 patients included in the analysis within 17 health facilities of the EsSalud Rebagliati Network
CKD Stages 3b-4
Overall
Outcome
Characteristic N = 2,798 Alive w/o Kidney Failure, N = 1,933 Death w/o Kidney Failure, N = 683 Kidney Failure, N = 182
Healthcare center
C. M. MALA 26 (100.0%) 10 (38.5%) 10 (38.5%) 6 (23.1%)
CAP II LURIN 12 (100.0%) 8 (66.7%) 4 (33.3%) 0 (0.0%)
CAP III LOS PROCEDES DE SJM 129 (100.0%) 87 (67.4%) 32 (24.8%) 10 (7.8%)
CAP III SAN ISIDRO 113 (100.0%) 68 (60.2%) 35 (31.0%) 10 (8.8%)
CAP III SURQUILLO 115 (100.0%) 81 (70.4%) 24 (20.9%) 10 (8.7%)
HOSP. II ANGAMOS 921 (100.0%) 644 (69.9%) 212 (23.0%) 65 (7.1%)
HOSPITAL EDGARDO REBAGLIATI MARTINS 22 (100.0%) 11 (50.0%) 8 (36.4%) 3 (13.6%)
HOSPITAL I CARLOS ALCANTARA B. 271 (100.0%) 194 (71.6%) 67 (24.7%) 10 (3.7%)
HOSPITAL I ULDARICO ROCCA F. 70 (100.0%) 47 (67.1%) 17 (24.3%) 6 (8.6%)
HOSPITAL II CANETE 90 (100.0%) 59 (65.6%) 28 (31.1%) 3 (3.3%)
POLICLINICO CHEQUEOS LARCO 15 (100.0%) 13 (86.7%) 1 (6.7%) 1 (6.7%)
POLICLINICO CHINCHA 408 (100.0%) 294 (72.1%) 93 (22.8%) 21 (5.1%)
POLICLINICO JUAN JOSE RODRIGUEZ LAZO 216 (100.0%) 134 (62.0%) 70 (32.4%) 12 (5.6%)
POLICLINICO PABLO BERMUDEZ 215 (100.0%) 152 (70.7%) 47 (21.9%) 16 (7.4%)
POLICLINICO PROCERES 144 (100.0%) 106 (73.6%) 30 (20.8%) 8 (5.6%)
POLICLINICO SANTA CRUZ 23 (100.0%) 18 (78.3%) 4 (17.4%) 1 (4.3%)
POLICLINICO VILLA MARIA 8 (100.0%) 7 (87.5%) 1 (12.5%) 0 (0.0%)
1 n (%)

6.8 Table S8

Mostrar código
# Selection of group of patients----
group <- c("Stages 3-4")

# Creacion de datsets para IA a 5 y 2 años
data %>% 
  filter(ckd_stage == group) -> data_filt

vdata.w <- crprep(
  Tstop = "time",
  status = "eventd",
  trans = c(1, 2),
  id = "id",
  keep = c("age", "male", "eGFR_mdrd", "acr", "risk2y", "risk5y"),
  data = data_filt
)

vdata.w1 <- vdata.w %>% filter(failcode == 1)
vdata.w2 <- vdata.w %>% filter(failcode == 2)

# For kidney failure
mfit_vdata1 <- survfit(
  Surv(Tstart, Tstop, status == 1) ~ 1,
  data = vdata.w1, 
  weights = weight.cens
)

smfit_vdata1 <- summary(mfit_vdata1, times = c(1, 2, 3, 4, 5))

res_ci_stg1 <- cbind(
  100 * (1 - smfit_vdata1$surv),
  100 * (1 - smfit_vdata1$upper),
  100 * (1 - smfit_vdata1$lower)
)

res_ci_stg1 <- round(res_ci_stg1, 2)

rownames(res_ci_stg1) <- c(
  "1-year", "2-year",
  "3-year", "4-year",
  "5-year"
)

colnames(res_ci_stg1) <- c(
  "Estimate", "Lower .95",
  "Upper .95"
)


# For death without kidney failure 
mfit_vdata2 <- survfit(
  Surv(Tstart, Tstop, status == 2) ~ 1,
  data = vdata.w2, 
  weights = weight.cens
)

smfit_vdata2 <- summary(mfit_vdata2, times = c(1, 2, 3, 4, 5))

res_ci_stg2 <- cbind(
  100 * (1 - smfit_vdata2$surv),
  100 * (1 - smfit_vdata2$upper),
  100 * (1 - smfit_vdata2$lower)
)

res_ci_stg2 <- round(res_ci_stg2, 2)

rownames(res_ci_stg2) <- c(
  "1-year", "2-year",
  "3-year", "4-year",
  "5-year"
)

colnames(res_ci_stg2) <- c(
  "Estimate", "Lower .95",
  "Upper .95"
)

# Selection of group of patients----
group <- c("Stages 3b-4")

# Creation of datasets for Cumulative Incidence at 2- and 5-yearsa
data %>% 
  filter(ckd_stage2 == group) -> data_filt

vdata.w <- crprep(
  Tstop = "time5y",
  status = "eventd5y",
  trans = c(1, 2),
  id = "id",
  keep = c("age", "male", "eGFR_mdrd", "acr", "risk2y", "risk5y"),
  data = data_filt
)

vdata.w1 <- vdata.w %>% filter(failcode == 1)
vdata.w2 <- vdata.w %>% filter(failcode == 2)

# For kidney failure
mfit_vdata3 <- survfit(
  Surv(Tstart, Tstop, status == 1) ~ 1,
  data = vdata.w1, 
  weights = weight.cens
)

smfit_vdata3 <- summary(mfit_vdata3, times = c(1, 2, 3, 4, 5))

res_ci_stg3 <- cbind(
  100 * (1 - smfit_vdata3$surv),
  100 * (1 - smfit_vdata3$upper),
  100 * (1 - smfit_vdata3$lower)
)

res_ci_stg3 <- round(res_ci_stg3, 2)

rownames(res_ci_stg3) <- c(
  "1-year", "2-year",
  "3-year", "4-year",
  "5-year"
)

colnames(res_ci_stg3) <- c(
  "Estimate", "Lower .95",
  "Upper .95"
)


# For death without kidney failure 
mfit_vdata4 <- survfit(
  Surv(Tstart, Tstop, status == 2) ~ 1,
  data = vdata.w2, 
  weights = weight.cens
)

smfit_vdata4 <- summary(mfit_vdata4, times = c(1, 2, 3, 4, 5))

res_ci_stg4 <- cbind(
  100 * (1 - smfit_vdata4$surv),
  100 * (1 - smfit_vdata4$upper),
  100 * (1 - smfit_vdata4$lower)
)

res_ci_stg4 <- round(res_ci_stg4, 2)

rownames(res_ci_stg4) <- c(
  "1-year", "2-year",
  "3-year", "4-year",
  "5-year"
)

colnames(res_ci_stg4) <- c(
  "Estimate", "Lower .95",
  "Upper .95"
)

# Table for 3a-3b-4 CKD patients----
res_ci <- cbind(res_ci_stg1, res_ci_stg2)

res_ci %>% 
  as_tibble(rownames = "Year") %>% 
  mutate(
    est1 = str_glue("{Estimate}%"), 
    est1.ci = str_glue(" ({`Lower .95`}% to {`Upper .95`}%)"), 
    est2 = str_glue("{V4}%"), 
    est2.ci = str_glue(" ({V5}% to {V6}%)")
  ) %>% 
  select(Year, est1, est1.ci, est2, est2.ci) %>% 
  mutate(across(everything(), ~as.character(.x))) -> res_ci2A 

res_ci2A %>% 
  flextable() %>% 
  add_header_row(values = c("Year", "Kidney failure", "Death without kidney failure"), colwidths = c(1, 2, 2)) %>% 
  set_header_labels(
    est1 = "%", 
    est1.ci = "95% CI", 
    est2 = "%", 
    est2.ci = "95% CI" 
  ) %>% 
  merge_v(j = 1, part = "header") %>% 
  set_caption("Table S8. Cumulative incidence of kidney failure and death without kidney failure in patients with CKD stages 3a-3b-4") %>% 
  theme_box() %>% 
  autofit() -> tab_cif1

tab_cif1 %>% 
  save_as_docx(path = here("Tables/TableS8.docx"))
Mostrar código
kable(
  res_ci2A, 
  caption = "Table S8. Cumulative incidence of Kidney Failure and Death w/o Kidney Failure in patients with CKD Stages 3a-3b-4", 
  col.names = c("", "%", "95%CI", "%", "95%CI")
  ) %>%
  kable_styling("striped", position = "center") %>%
  add_header_above(c(" " = 1, "Kidney Failure" = 2, "Death w/o Kidney Failure" = 2))
Table S8. Cumulative incidence of Kidney Failure and Death w/o Kidney Failure in patients with CKD Stages 3a-3b-4
Kidney Failure
Death w/o Kidney Failure
% 95%CI % 95%CI
1-year 0.68% (0.49% to 0.86%) 3.82% (3.38% to 4.25%)
2-year 1.52% (1.24% to 1.79%) 7.49% (6.89% to 8.08%)
3-year 2.23% (1.9% to 2.57%) 11.2% (10.48% to 11.91%)
4-year 2.88% (2.5% to 3.26%) 15.58% (14.75% to 16.41%)
5-year 3.37% (2.95% to 3.8%) 20.46% (19.48% to 21.42%)

The estimated cumulative incidence of renal failure at 2 and 5 years was 1.52% ( (1.24% to 1.79%)) and 3.37% ( (2.95% to 3.8%)) in patients with stage 3a, 3b or 4 CKD, respectively (Table S8).

6.9 Tabla S9

Mostrar código
# Tabla for 3b-4----
res_ci <- cbind(res_ci_stg3, res_ci_stg4)

res_ci %>% 
  as_tibble(rownames = "Year") %>% 
  mutate(
    est1 = str_glue("{Estimate}%"), 
    est1.ci = str_glue(" ({`Lower .95`}% to {`Upper .95`}%)"), 
    est2 = str_glue("{V4}%"), 
    est2.ci = str_glue(" ({V5}% to {V6}%)")
  ) %>% 
  select(Year, est1, est1.ci, est2, est2.ci) %>% 
  mutate(across(everything(), ~as.character(.x))) -> res_ci2B 

res_ci2B %>% 
  flextable() %>% 
  add_header_row(values = c("Year", "Kidney failure", "Death without kidney failure"), colwidths = c(1, 2, 2)) %>% 
  set_header_labels(
    est1 = "%", 
    est1.ci = "95% CI", 
    est2 = "%", 
    est2.ci = "95% CI" 
  ) %>% 
  merge_v(j = 1, part = "header") %>% 
  set_caption("Table S9. Cumulative incidence of kidney failure and death without kidney failure in patients with CKD stages 3b-4") %>% 
  theme_box() %>% 
  autofit() -> tab_cif2

tab_cif2 %>% 
  save_as_docx(path = here("Tables/TableS9.docx"))
Mostrar código
kable(
  res_ci2B, 
  caption = "Table S9. Cumulative incidence of Kidney Failure and Death w/o Kidney Failure in patients with CKD Stages 3b-4", 
  col.names = c("", "%", "95%CI", "%", "95%CI")
  ) %>%
  kable_styling("striped", position = "center") %>%
  add_header_above(c(" " = 1, "Kidney Failure" = 2, "Death w/o Kidney Failure" = 2))
Table S9. Cumulative incidence of Kidney Failure and Death w/o Kidney Failure in patients with CKD Stages 3b-4
Kidney Failure
Death w/o Kidney Failure
% 95%CI % 95%CI
1-year 1.54% (1.08% to 1.99%) 6% (5.12% to 6.88%)
2-year 3.15% (2.5% to 3.79%) 10.72% (9.57% to 11.86%)
3-year 4.72% (3.93% to 5.5%) 15.51% (14.16% to 16.84%)
4-year 6.02% (5.12% to 6.9%) 21.03% (19.48% to 22.55%)
5-year 6.86% (5.89% to 7.83%) 26.59% (24.84% to 28.3%)

In patients with stage 3b-4 CKD, the cumulative incidences at 2 and 5 years were 10.72% ( (9.57% to 11.86%)) and 26.59% ( (24.84% to 28.3%)), respectively. Table S9 details the cumulative incidences of renal failure from 1 to 5 years in both cohorts.

6.10 Fig S1

Mostrar código
data %>%
  filter(ckd_stage == "Stages 3-4") %>% 
  mutate(evento = "CKD 3-4") -> dataA

data %>%
  filter(ckd_stage2 == "Stages 3b-4") %>% 
  mutate(evento = "CKD 3b-4") -> dataB  

dataA %>% 
  bind_rows(dataB) -> data2
  
data2 %>%
  ggplot(aes(x = evento, y = age, fill = evento)) + 
  sm_raincloud() + 
  labs(x = "", 
       y = "Age, years") + 
  coord_cartesian(ylim = c(0, 100)) -> p1

data2 %>%
  ggplot(aes(x = evento, y = eGFR_ckdepi, fill = evento)) + 
  sm_raincloud() + 
  labs(x = "", y = "eGFR, ml/min/1.73m<sup>2</sup>") +  
  coord_cartesian(ylim = c(0, 60)) + 
  theme(axis.title.y = element_markdown()) -> p2

data2 %>%
  ggplot(aes(x = evento, y = acr, fill = evento)) + 
  sm_raincloud() + 
  labs(x = "", 
       y = "ACR, mg/g") + 
  coord_cartesian(ylim = c(0, 150000)) -> p3

data2 %>%
  ggplot(aes(x = evento, y = acr, fill = evento)) + 
  sm_raincloud() + 
  scale_y_continuous(trans = "log10", labels = label_log()) +
  labs(x = "", 
       y = "ACR, mg/g (log scale)") -> p4

(p1 + p2 ) / (p3 + p4 ) + plot_annotation(tag_levels = "A") -> plot_dist_vars

ggsave(filename = "plot_dist_vars.png", 
       plot = plot_dist_vars, 
       device = "png",
       path = here("Figures/"), 
       scale = 2, 
       width = 18, 
       height = 18, 
       units = "cm",
       dpi = 900
       )
Mostrar código
knitr::include_graphics(here("Figures/plot_dist_vars.png"))

Figure 3: Distribution of four variables of KFRE equation in CKD 3-4 and CKD 3b-4 patientes. (A) age in years, (B) estimate glomerural filtration rate (eGFR) according to CKD-EPI formula, (C) urine albumin to creatinine ratio (ACR) expressed in original scale and in (D) natural logarithm scale for better comparison of the distributions.

6.11 Fig S2

Mostrar código
data %>%
  filter(ckd_stage == "Stages 3-4") %>% 
  mutate(
    evento2y = case_when(
      eventd2ylab %in% c("Alive w/o Kidney Failure", 
                         "Death w/o Kidney Failure") ~ "No", 
      eventd2ylab == "Kidney Failure" ~ "Yes", 
      TRUE ~ as.character(NA)), 
    evento5y = case_when(
      eventd5ylab %in% c("Alive w/o Kidney Failure", 
                         "Death w/o Kidney Failure") ~ "No", 
      eventd5ylab == "Kidney Failure" ~ "Yes", 
      TRUE ~ as.character(NA))
    ) -> dataA

data %>%
  filter(ckd_stage2 == "Stages 3b-4") %>% 
  mutate(
    evento2y = case_when(
      eventd2ylab %in% c("Alive w/o Kidney Failure", 
                         "Death w/o Kidney Failure") ~ "No", 
      eventd2ylab == "Kidney Failure" ~ "Yes", 
      TRUE ~ as.character(NA)), 
    evento5y = case_when(
      eventd5ylab %in% c("Alive w/o Kidney Failure", 
                         "Death w/o Kidney Failure") ~ "No", 
      eventd5ylab == "Kidney Failure" ~ "Yes", 
      TRUE ~ as.character(NA))
    ) -> dataB  
  
dataA %>%
  ggplot(aes(x = risk2y, color = evento2y, fill = evento2y)) + 
  geom_histogram(alpha = 0.4) + 
  labs(x = "2-years predicted risk by KFRE", 
       y = "Frequency count", 
       color = "Kidney failure", 
       fill = "Kidney failure", 
       title = "CKD stages 3-4") + 
  theme_classic() + 
  theme(plot.title = element_text(hjust = 0.5, 
                                  face = "bold")) + 
  coord_cartesian(xlim = c(0, 1)) -> p1

dataA %>%
  ggplot(aes(x = risk5y, color = evento5y, fill = evento5y)) + 
  geom_histogram(alpha = 0.4) + 
  labs(x = "5-years predicted risk by KFRE", 
       y = "Frequency count", 
       color = "Kidney failure", 
       fill = "Kidney failure") + 
  theme_classic() + 
  coord_cartesian(xlim = c(0, 1)) -> p2

dataB %>%
  ggplot(aes(x = risk2y, color = evento2y, fill = evento2y)) + 
  geom_histogram(alpha = 0.4) + 
  labs(x = "2-years predicted risk by KFRE", 
       y = "Frequency count", 
       color = "Kidney failure", 
       fill = "Kidney failure", 
       title = "CKD stages 3b-4") + 
  theme_classic() + 
  theme(plot.title = element_text(hjust = 0.5, 
                                  face = "bold")) + 
  coord_cartesian(xlim = c(0, 1)) -> p3

dataB %>%
  ggplot(aes(x = risk5y, color = evento5y, fill = evento5y)) + 
  geom_histogram(alpha = 0.4) + 
  labs(x = "5-years predicted risk by KFRE", 
       y = "Frequency count", 
       color = "Kidney failure", 
       fill = "Kidney failure") + 
  theme_classic() + 
  coord_cartesian(xlim = c(0, 1)) -> p4

((p1 / p2) | (p3 / p4)) +  
  plot_layout(guides = "collect") + 
  plot_annotation(tag_levels = "A")  -> plot_dist_risks

ggsave(filename = "plot_dist_risks.png", 
       plot = plot_dist_risks, 
       device = "png",
       path = here("Figures/"), 
       scale = 2, 
       width = 18, 
       height = 18, 
       units = "cm",
       dpi = 900
       )
Mostrar código
knitr::include_graphics(here("Figures/plot_dist_risks.png"))

Figure 4: Distribution of the 2-year and 5-year predicted risk estimated by KFRE equation according kidney failure

6.12 Table S13

Mostrar código
# Soure: Table 1, subtotal of Non-North America (doi: 10.1001/jama.2015.18202)
table_initial_kfre <- data.frame(
  vars = c("Numbef of participants", 
           "F/U Time, years, Median (IQR)", 
           "Age, years (SD)", 
           "Male, n (%)", 
           "Black ethnicity, n (%)", 
           "eGFR, ml/min/1.73m2 (SD)", 
           "Albuminuria, n (%)", 
           "Kidney Failure Incidence (per 1000 py)"), 
  original = c("103753", 
               "4 (3, 6)", 
               "71 (12)", 
               "46632 (45%)", 
               "393 (0.4%)",  
               "47 (12)", 
               "24962 (34%)", 
               "9.2"), 
  ckdA = c("7519", 
           "4.9 (3.5, 5.9)", 
           "74 (10.2)", 
           "3412 (45.4%)", 
           "0 (0%)", 
           "46.2 (9.8)", 
           "2747 (36.5%)", 
           "7.4"), 
  ckdB = c("2798", 
           "4.6 (3.2, 5.8)", 
           "75.6 (10.6)", 
           "1400 (50%)", 
           "0 (0%)", 
           "35.7 (7.3)", 
           "1304 (46.6%)", 
           "16.1")
  ) 

table_initial_kfre %>% 
  gt() %>% 
  cols_label(vars = "Characteristics", 
             original = "Original Study (Non-North American population)", 
             ckdA = "CKD Stages 3a-4", 
             ckdB = "CKD Stages 3b-4") %>% 
  tab_spanner(label = "Current study", 
              columns = c(ckdA, ckdB)) -> tabS13

tabS13 %>% 
  gtsave(filename = "TableS13.docx", path = here("Tables"))

tabS13
Characteristics Original Study (Non-North American population) Current study
CKD Stages 3a-4 CKD Stages 3b-4
Numbef of participants 103753 7519 2798
F/U Time, years, Median (IQR) 4 (3, 6) 4.9 (3.5, 5.9) 4.6 (3.2, 5.8)
Age, years (SD) 71 (12) 74 (10.2) 75.6 (10.6)
Male, n (%) 46632 (45%) 3412 (45.4%) 1400 (50%)
Black ethnicity, n (%) 393 (0.4%) 0 (0%) 0 (0%)
eGFR, ml/min/1.73m2 (SD) 47 (12) 46.2 (9.8) 35.7 (7.3)
Albuminuria, n (%) 24962 (34%) 2747 (36.5%) 1304 (46.6%)
Kidney Failure Incidence (per 1000 py) 9.2 7.4 16.1

7 Ticket de Reprocubilidad

Mostrar código
R version 4.2.3 (2023-03-15 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Peru.utf8  LC_CTYPE=Spanish_Peru.utf8   
[3] LC_MONETARY=Spanish_Peru.utf8 LC_NUMERIC=C                 
[5] LC_TIME=Spanish_Peru.utf8    

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] smplot2_0.1.0             janitor_2.2.0            
 [3] gt_0.8.0                  DiagrammeR_1.0.9         
 [5] ggtext_0.1.2              htmltools_0.5.4          
 [7] devtools_2.4.5            usethis_2.1.6            
 [9] gghalves_0.1.4            xml2_1.3.3               
[11] downlit_0.4.2             broom_1.0.3              
[13] dcurves_0.4.0             glue_1.6.2               
[15] labelled_2.10.0           scales_1.2.1             
[17] cowplot_1.1.1             ggsci_2.9                
[19] survminer_0.4.9           ggpubr_0.6.0             
[21] patchwork_1.1.2           webshot_0.5.4            
[23] gridExtra_2.3             rsample_1.1.1            
[25] lubridate_1.9.2           forcats_1.0.0            
[27] stringr_1.5.0             dplyr_1.1.0              
[29] purrr_1.0.1               readr_2.1.4              
[31] tidyr_1.3.0               tibble_3.1.8             
[33] tidyverse_2.0.0           boot_1.3-28.1            
[35] gtsummary_1.7.0           flextable_0.8.6          
[37] kableExtra_1.3.4          knitr_1.42               
[39] plotrix_3.8-2             pec_2022.05.04           
[41] prodlim_2019.11.13        pseudo_1.4.3             
[43] geepack_1.3.9             KMsurv_0.1-5             
[45] mstate_0.3.2              riskRegression_2022.11.28
[47] cmprsk_2.2-11             rms_6.5-0                
[49] SparseM_1.81              Hmisc_4.8-0              
[51] ggplot2_3.4.1             Formula_1.2-5            
[53] lattice_0.20-45           survival_3.5-3           
[55] skimr_2.1.5               here_1.0.1               

loaded via a namespace (and not attached):
  [1] pacman_0.5.1            utf8_1.2.3              tidyselect_1.2.0       
  [4] htmlwidgets_1.6.1       grid_4.2.3              munsell_0.5.0          
  [7] codetools_0.2-19        ragg_1.2.5              interp_1.1-3           
 [10] miniUI_0.1.1.1          future_1.31.0           withr_2.5.0            
 [13] colorspace_2.1-0        highr_0.10              uuid_1.1-0             
 [16] rstudioapi_0.14         ggsignif_0.6.4          officer_0.6.0          
 [19] fontLiberation_0.1.0    listenv_0.9.0           labeling_0.4.2         
 [22] repr_1.1.6              mets_1.3.2              farver_2.1.1           
 [25] rprojroot_2.0.3         parallelly_1.34.0       vctrs_0.5.2            
 [28] generics_0.1.3          TH.data_1.1-1           xfun_0.37              
 [31] timechange_0.2.0        fontquiver_0.2.1        markdown_1.5           
 [34] R6_2.5.1                timereg_2.0.5           cachem_1.0.7           
 [37] promises_1.2.0.1        multcomp_1.4-22         nnet_7.3-18            
 [40] gtable_0.3.1            globals_0.16.2          processx_3.8.0         
 [43] sandwich_3.0-2          rlang_1.0.6             MatrixModels_0.5-1     
 [46] systemfonts_1.0.4       rstatix_0.7.2           checkmate_2.1.0        
 [49] yaml_2.3.7              abind_1.4-5             backports_1.4.1        
 [52] httpuv_1.6.9            gridtext_0.1.5          tools_4.2.3            
 [55] lava_1.7.2.1            ellipsis_0.3.2          RColorBrewer_1.1-3     
 [58] sessioninfo_1.2.2       Rcpp_1.0.10             visNetwork_2.1.2       
 [61] base64enc_0.1-3         prettyunits_1.1.1       ps_1.7.2               
 [64] rpart_4.1.19            openssl_2.0.5           deldir_1.0-6           
 [67] urlchecker_1.0.1        zoo_1.8-11              haven_2.5.2            
 [70] cluster_2.1.4           fs_1.6.1                furrr_0.3.1            
 [73] crul_1.3                magrittr_2.0.3          data.table_1.14.8      
 [76] mvtnorm_1.1-3           pkgload_1.3.2           hms_1.1.2              
 [79] mime_0.12               evaluate_0.20           xtable_1.8-4           
 [82] jpeg_0.1-10             compiler_4.2.3          fontBitstreamVera_0.1.1
 [85] crayon_1.5.2            later_1.3.0             tzdb_0.3.0             
 [88] MASS_7.3-58.2           broom.helpers_1.12.0    sdamr_0.2.0            
 [91] Matrix_1.5-3            car_3.1-1               cli_3.6.0              
 [94] parallel_4.2.3          pkgconfig_2.0.3         km.ci_0.5-6            
 [97] numDeriv_2016.8-1.1     foreign_0.8-84          foreach_1.5.2          
[100] svglite_2.1.1           rvest_1.0.3             snakecase_0.11.0       
[103] callr_3.7.3             digest_0.6.31           httpcode_0.3.0         
[106] rmarkdown_2.20          survMisc_0.5.6          htmlTable_2.4.1        
[109] gdtools_0.3.1           curl_5.0.0              commonmark_1.8.1       
[112] shiny_1.7.4             quantreg_5.94           pwr_1.3-0              
[115] lifecycle_1.0.3         nlme_3.1-162            jsonlite_1.8.4         
[118] carData_3.0-5           viridisLite_0.4.1       askpass_1.1            
[121] fansi_1.0.4             pillar_1.8.1            pkgbuild_1.4.0         
[124] fastmap_1.1.1           httr_1.4.5              remotes_2.4.2          
[127] zip_2.2.2               png_0.1-8               iterators_1.0.14       
[130] sass_0.4.5              profvis_0.3.7           stringi_1.7.12         
[133] polspline_1.1.22        textshaping_0.3.6       gfonts_0.2.0           
[136] latticeExtra_0.6-30     memoise_2.0.1           future.apply_1.10.0